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ABSTRACT 


A three-dimensional computer code (KIVA) is being modified to include state-of- 
the-art submodels for diesel engine flow and combustion: spray atomization, drop 
breakup/coalescence, multi-component fuel vaporization, spray /wall interaction, 
ignition and combustion, wall heat transfer, unburned HC and NOx formation, 
soot and radiation and the intake flow process. 

Improved and/or new submodels which have been completed are: wall heat 
transfer with unsteadiness and compressibility, laminar-turbulent characteristic 
time combustion with unburned HC and Zeldo’vich NOx, and spray/ wall 
impingement with rebounding and sliding drops. Results to date show that 
adding the effects of unsteadiness and compressibility improves the accuracy of 
heat transfer predictions; spray drop rebound can occur from walls at low 
impingement velocities (e.g., in cold-starting); larger spray drops are formed at the 
nozzle due to the influence of vaporization on the atomization process; a laminar- 
and-turbulent characteristic time combustion model has the flexibility to match 
measured engine combustion data over a wide range of operating conditions; and, 
finally, the characteristic time combustion model can also be extended to allow 
predictions of ignition. 


The accuracy of the predictions is being assessed by comparisons with available 
measurements. Additional supporting experiments are also described briefly. To 
date, comparisons have been made with measured engine cylinder pressure and 
heat flux data for homogeneous charge, spark-ignited and compression-ignited 
engines, and also limited comparisons for diesel engines. The model results are in 
good agreement with the experiments. 
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INTRODUCTION 


A detailed understanding of diesel engine combustion is required in order to 
work effectively at improving performance and reducing emissions while not 
compromising fuel economy. The objective of this research is to apply 
advanced modeling techniques to study the influence of in-cylinder processes 
on efficiency and pollutant emissions. The three-dimensional code, KTVA 
[1,2], has been selected for use since it is the most developed of available codes. 

State-of-the-art submodels for the important physical processes in diesel 
combustion are being included in KIVA as part of the research effort. This 
report summarizes progress to date on models for wall heat transfer, fuel 
drop vaporization, spray vaporization, atomization, ignition and the intake 
flow process. Future activities under the grant are also described briefly . 

Close contact is being maintained with engine industries during the course of 
the research. This includes participation in DOE diesel engine working group 
meetings 1 . In addition, technology has been transferred directly to industry 
(e.g., we have already made the FORTRAN subroutines listed in Appendix 3 
of this proposal available to the Software Development Group, Engine 
Division, Caterpillar Inc.). 

Other activities that have helped the research effort include our pro-active 
role in the formation and organization of the KIVA users group (which 
currently numbers about 80 organizations in the U.S. - see Appendix 1). This 
activity has had the additional benefit of keeping us informed of code 
enhancements and new developments elsewhere in the user community. 
Improvements are being made to the KIVA code itself at the Los Alamos 
National Laboratories. We have been able to take advantage of their 
enhancements by ensuring that our submodels are transportable in the form 
of subroutines (e.g.. Appendix 3). 

BACKGROUND 

The diesel engine is the leading heavy-duty power plant because of its 
superior energy efficiency. However, because of environmental concerns, 
proposed federal emission standards require reductions in both nitric oxides 
(NOx) and particulates for heavy-duty diesel engines. A detailed 
understanding of combustion is required in order to work effectively at 
reducing these by-products of combustion within the engine cylinder, while 
still not compromising engine fuel economy. Alternative methods of 


Prof. Reitz gave the presentations 'Modeling Diesel Engine Spray Evaporation arid 
Combustion' at the Spring 1991 DOE Diesel Group Meeting, Cummins Engine Company, 
Indiana, May 2, 1991, and '3-D Modeling of Diesel Spray Vaporization and Combustion,' 
Fall 1991 DOE Diesel Group Meeting, Pennsylvania State University, October 23, 1991. 
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controlling emissions outside of the engine cylinder by means of particulate 
traps and catalysts are not appealing due to their high complexity and cost. 

Electronically controlled diesel fuel injection equipment has been under 
active development for the past few years. Recent developments in 
microprocessor technology have created opportunities for better control of 
fuel delivery and timing than has hitherto been possible. This improved 
flexibility over the injection parameters offers the possibility of major 
improvements in diesel engine performance [3]. Improved injection timing 
flexibility and fast closing times are factors that are known to impact engine 
emissions [4,5]. However, the improved flexibility also introduces more 
configuration possibilities, and this complicates the task of selecting optimum 
injection parameters. 

The objective of this research program is to apply advanced modeling 
techniques to provide guidance for the selection of fuel injection and 
combustion parameters. This will help in the evolution of more fully 
optimized, clean and efficient diesel engines. Three-dimensional 
mathematical models such as the KIVA code are now available and offer a 
means to help gain the needed understanding of engine combustion and can 
accelerate the rate of progress in engine development beyond that obtained 
with traditional 'cut and try' approaches. 


In the compression ignition (Cl) diesel engine a number of steps are involved 
in the combustion process. The air is raised to a high temperature on the 
compression stroke. Multiple high pressure jets of fuel are injected into the 
chamber which disintegrate into dense vaporizing sprays. After some delay 
time, self-ignition occurs at points within the sprays. A period of rapid 
pressure rise follows, during which the vaporized (premixed) portions of the 
fuel are burned. In the final stages of combustion, diffusion and mixing 
processes lead to a controlled rate of pressure rise. In some engines, part of 
the injected fuel impinges on the combustion chamber walls and the final 
stage of combustion is then influenced by the rate of vaporization of the fuel 
from the walls. Thus, engine performance (both efficiency and emissions 
levels) depends on the spatial and temporal details of the atomization of the 
liquid, vaporization of the droplets, mixing of the fuel and air, and other 
factors. It is necessary to understand all of these processes and to identify their 
controlling parameters in order to improve engine performance. 

There is evidence that in-cylinder control of pollutant emission levels is 
possible through the use of ultra-high injection pressures and modified 
injector nozzle designs or modified in-cylinder gas flow-fields [6,7]. Also, 
optimized low-heat-rejection (LHR) engines have demonstrated significant 
reductions in smoke and particulates [8]. Unfortunately, the high in-cylinder 
gas temperatures of LHR engines also lead to high NOx levels and methods 
for NOx control such as injection timing retard tend to increase fuel 
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consumption. It has been difficult to achieve satisfactory pollutant control 
strategies over wide ranges of engine operation with these approaches alone. 

The present advanced modeling techniques will help supply engine and fuel 
injection design guidelines in this study. The approach taken is to use a 
three-dimensional code that includes the most advanced available submodels 
for the important processes that influence diesel combustion. The model 
accounts for the complex interactions that occur in the turbulent, multiphase, 
combusting three-dimensional flow in the engine. The complexity of these 
interactions has made it difficult to gain the necessary level of understanding 
required for in-cylinder control of emissions by traditional methods. 

RESEARCH PROGRAM 

The research program was initiated in the Fall of 1989 as a five year program. 
The emphasis of the research is on the application of a comprehensive engine 
combustion code to assess the effect of the interacting subprocesses on diesel 
engine performance, rather than on the development of new models for the 
subprocesses themselves. The elements of the code are being assembled from 
existing state-of-the-art submodels. This use of multidimensional modeling 
as an engine development tool is timely and justifiable due to recent 
advances in submodel formulations. 

The program was expanded to include the modeling of the diesel engine 
intake process in 1990 due to its importance in the combustion process. A 
summary of the 1989-1990 research activities has been described in Refs. [12- 
14]. The research has been organized into three main activities: Part A - 
Submodel Implementation, Part B - Model Application and Part C - 
Experiments. Part A tasks are described in the Progress-to-date section 
(page 7). Part B tasks were started in the the summer of 1991 and are described 
on page 33. The Part C tasks are described briefly in the Future Work section 
of this report on page 35. 

KIVA CODE 

The KIVA code has been chosen for use in the present work since it is the 
most developed of the available 3-D engine models. The code solves the 
conservation equations for the transient dynamics of vaporizing fuel sprays 
interacting with flowing multicomponent gases which are undergoing 
mixing, ignition, chemical reactions, and heat transfer. The code has the 
ability to calculate flows in engine cylinders with arbitrary shaped piston 
geometries, including the effects of turbulence and wall heat transfer. 

Accuracy of KIVA predictions has been assessed by comparisons with 
measured cylinder pressures in direct-injection diesel [9], direct-injection 
stratified-oharge rotary engines [10], and in homogeneous-charge reciprocating 
engines [11]. 
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Due to practical limitations of computer storage and run times it is necessary 
to introduce submodels into engine codes for processes that occur on time 
and length scales that are too short to be resolved (e.g., to resolve the flow- 
field around 10 |im diameter drops in a 10 cm diameter combustion chamber 
requires about 10^ grid points. A practical limit for current computers is 
about 10 5 grid points). The use of submodels to describe unresolved physical 
processes introduces some empiricism into the computations. However, the 
compromise between accuracy and feasibility of computation is justified by 
the insight which model calculations offer and confidence in the model 
predictions is gained by comparison with experiments. 

In the present work, the KIVA code is being modified to include existing 
state-of-the-art submodels for important processes occurring in diesel engines. 
Some of the submodels already appear in the standard KIVA-II code (e.g., the 
drop drag, turbulence dispersion and turbulence modulation submodels - see 
How Chart of KIVA Subroutines below). Others have already been 
developed and applied in earlier separate versions of KIVA. Thus, the 
present research also serves to assemble the most recent models into one 
comprehensive version of KIVA. The Subroutines that have been modified 
to data in this effort are identified with asterisks in the flow chart. 

PROGRESS TO-DATE 

Progress to date in each of the Part A tasks scheduled for the first and second 
years is discussed in this section and summarized in Table 1 (page 9). 
Additional details of progress are given in Ref. [14]. Appendix 2 gives 
additional details of published papers connected with the research activity. 
Appendix 3 contains FORTRAN source code listings for the spray/ wall 
impingement, combustion and atomization submodels that have been 
upgraded as part of the research activity. 

1. Implement KIVA on TITAN 

Due to the high cost of super-computer time, it is important to demonstrate 
that the KIVA code can be run on dedicated super- workstations. This will 
ensure that the code is attractive to industrial users. The KIVA-2 code has 
been running on the Stardent TITAN computer at the Engine Research 
Center since 9/89. With the model P-3 CPU boards (installed 8/90) we are able 
to vectorize KIVA-2 satisfactorily. Our benchmark runs indicate that the 
TITAN now runs only a factor of 10 slower than a CRAY XMP. Also, under 
ARO funding, the TITAN has been expanded from two to three CPUs, which 
has increased productivity. In addition, we have been able to secure 
donations of CRAY computer time from CRAY Research, Inc. and from the 
San Diego Supercomputer Installation. Thus, the availability of computer 
time has not impeded progress on the grant. 
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FLOW CHART of KIVA SUBROUTINES Refs. [1,2] 

Routines to be upgraded in the present effort are indicated by the asterisks. 
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A. Submodel Implementation 


Table 1 - Submodel Implementation 


Year 

Phase / Activity 

9/1/89 

9/1/90 

1 

9/1/90 

9/1/91 

2 

9/1/91 

9/1/92 

3 


1*14 Wi*£ 

1 

Implement KIVA-2 on 
TITAN 

X 

X 




2 

Atomization, 
spray/wall interaction, 
heat transfer 

X 





3 

HC and NOx emission 

X 





4 

Multicomponent fuel 
Vaporization 

X 

X 




5 

Diesel ignition and 
combustion 


X 

X 



6 

Soot and radiation 



X 

X 


7 

Intake flow modeling 
Phase 1 Mesh 
Phase 2 Mesh 
Moving valves 


X 

X 

X 
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2a. Spray/Wall Interaction 

The KIVA-2 code does not contain a spray/ wall interaction model and the 
model of Naber and Reitz [15] has been included in KIVA. This model has 
been refined further to account for rebounding and sliding drops [16] (see 
Appendix 3). Upon impact, low Weber number drops rebound from the 
surface with an angle of rebound that is determined by curve-fitting 
experimental data on single drops. The need to consider drop rebound was 
discovered when computations of diesel engine cold-starting were made [16]. 

In this case, drop rebound provides a mechanism for fuel (vapor) to penetrate 
back into the central regions of the chamber where the gas temperature is 
high enough for ignition to occur. Without drop rebound the fuel is 
predicted to accumulate near the walls where the gas temperature is too low 
for ignition [16]. 

2b. Wall Heat Transfer in Premixed-Charge Engine Combustion 

Wall heat transfer models have been tested for premixed-charge, spark-ignited- 
engine combustion as described in Ref [17]. Since the combustion model in 
KIVA does not account for the influence of turbulence on combustion, a laminar 
and turbulent combustion submodel has been added for this study (see also task 5 
below). 

Typical results showing the application of the model to spark-ignited premixed- 
charge engine combustion are shown in Fig. 1 [17]. As can be seen in Fig. la, the 
computed results (squares) reproduce the form of measured pressure data 
(circles) very well. 



Crank angle (degrees) 



-30 -20 -10 0 10 20 30 

Crank angle (degrees) 


Figure 1 Comparison between measured and predicted cylinder pressure 
a.) - left, and wall heat flux b.) - right, for premixed-charge engine combustion 
using the characteristic time combustion model and a heat transfer model that 
includes the effect of compressibility and unsteadiness (modified case). 

Engine speed is 1500 rev/min. 




1 0 





The wall heat flux predictions in Fig. lb were obtained using the original 
KIVA 'log-law-of-the-wall' heat transfer model (squares) and a modified 
model that accounts for the effects of compressibility and unsteadiness on the 
heat transfer (triangles) [17]. The predicted heat transfer is seen to be much 
improved using the upgraded model, indicating that unsteadiness and 
compressibility effects need to be considered in engine computations. 
However it is also evident that additional improvements are needed in heat 
transfer models, since there is still a 10-20% discrepancy between measured 
and predicted heat flux values [17]. Reference [17] describes other results that 
test the wall heat transfer and combustion models over a range of engine 
operating conditions (motored and fired conditions at two engine speeds). 

Further work is in progress on testing the heat transfer models using 
homogeneous charge, compression ignition engine data (see Task 5, below). 

2c. Atomization, vaporization and combustion 

Spray computations have been made using our atomization model that is 
based on considerations of surface wave instabilities [18]. In modeling the 
atomization process in high pressure sprays, it has been shown that the drop 
size and penetration length predicted by KIVA are in good agreement with 
experimental data for non-evaporating sprays at room temperature [18]. 

However, for high pressure and high temperature sprays, significant 
discrepancies between the computational and experimental data have been 
observed [19]. As shown in Fig. 2, using the grid resolution that ensures grid 
free results in non-vaporizing sprays, the computed vapor penetrations 
under-estimate the measurements. Since spray penetration directly 
influences combustion [19], heat transfer, soot formation and other processes, 
identifying the causes of this discrepancy and making corresponding 
modifications in KIVA has become essential to continued progress on the 
grant. 

The computations in Fig. 2 were compared with experimental data of 
Kamimoto et al. [20] for room temperature liquid n-tridecane (C13H28) 
injections from a single hole nozzle (d=0.16mm) into nitrogen gas at 850 - 900 
K. The gas pressure was maintained at 3 MPa while the injection pressure 
were varied from 30 MPa to 110 MPa. 

The effects of several factors were considered in the present study, including 
the influence of vaporization and gas compressibility on the atomization 
process, and the effects of model constants and numerical parameters. 

The spray penetration was found to be very sensitive to the grid size, 
especially to the grid size near the exit of nozzle, as indicated in Fig. 3. For 
relatively coarse meshes the penetration increases gradually as the grid size 
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Figure 2 Computed spray vapor penetration versus time. Symbols are 
measurements of Kamimoto [20] and lines are KTVA predictions (penetration 
is defined as that distance from the nozzle exit where the fuel vapor 
concentration reaches 10% of the maximum value in the spray). 



Figure 3 Effect of numerical mesh grid size on spray penetration. AP=30 MPa, 
Line 1: grid size=l x 5 mm 

Line 2: non-uniform grid 1.25 x 2 mm (near nozzle) to 1.25 x 4.5mm 
Line 3: uniform grid 1.25 x 2 mm. 

Line 4: non-uniform grid 0.75 x 1 mm (near nozzle) to 0.75 x 3 mm 
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decreases, which explains why the effect was not noticed in our previous 
studies [19]. For very fine meshes the penetration is increased dramatically. 
An explanation for the observed mesh-dependence is that the fuel vaporized 
from a drop is distributed evenly over the whole computational cell that 
contains that drop at the end of each timestep in KIVA. This procedure 
might be reasonable for very fine meshes (or large timesteps) but is hardly 
satisfied for coarse meshes where the vapor would not have enough time to 
diffuse or mix over the whole cell. Hence for coarse meshes the local vapor 
concentration around a drop should be higher than the cell-volume-averaged 
value. This leads to an over-estimate of the drop evaporation rate and a 
corresponding under-estimate in the spray penetration. 

The implications of this sensitivity to grid resolution for engine 
computations, where computer cost considerations limit grid sizes, are 
substantial. One solution to the grid sensitivity problem could be to develop 
correction procedures for use with coarse grids. For example, assume that 
instead of occupying the whole cell, the vapor evaporated from a drop is 
confined within a sphere with radius ~ (D t) 1/2 , where D is the vapor 
diffusivity in the gas and t is the lifetime of the drop. The vapor mass 
fraction is now averaged over the sphere instead of over the entire cell (until 
the sphere volume reaches that of the cell, see Fig. 4). Figure 5 shows the 
penetration thus calculated for a coarse mesh (1.25 x 2 mm) which can be seen 
to be much improved. 


Sub-grid-scale vapor diffusion model 



Figure 4 Schematic diagram of the sub-grid-scale vapor diffusion model. The 
fuel vapor concentration is calculated using the volume of the sphere with 
radius R instead of the cell volume until the vapor fills the cell. 
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Line 1: cell- volume averaged vapor (standard KIVA with - 1.25 x 2 mm mesh) 
Line 2 : new vapor-sphere averaged correction method ( 1.25 x 2 mm mesh) 
Line 3 : same as line 1, but with 0.75 x 1 mm mesh. 

The comparisons with the constant volume bomb experiments demonstrated 
the importance of grid resolution in vaporizing spray calculations. For that 
reason, the engine calculations have been made using fine grids in an attempt 
to accurately resolve the spray. The computations were made for Cummins 
NH engine for which extensive experimental data are available with 
measurements of injection characteristics, cylinder pressure, and flame 
temperatures from Yan and Borman [21]. Results using a cycle analysis 
simulation program are available for the same engine [22]. Data, such as the 
wall temperatures, calculated using the steady state mode were used as input 
for the initial conditions in the multidimensional calculations as shown in 
the table below. Additional details of the computations are described in 
Gonzalez et al. [23]. 

The calculations were started at inlet valve closing (150 deg BTDC) and 
considered a 45 degree sector of the engine that included one of the eight 
spray plumes (i.e., sector symmetry was assumed). The results shown in Figs. 

6 and 7 compare spray penetration computations (made without considering 
combustion) for two different grids. The first grid used 25x6x18 cells with 
equal spacing in the azimuthal direction. The second grid used the same 
number of axial cells but had nonuniform azimuthal spacing and gives much 
increased resolution in the vicinity of the spray - the minimum azimuthal 
angle is 2 degrees compared to 7.5 degrees in the equally spaced grid. Both 
grids had 7 axial planes in the bowl and 3 planes in the squish area at TDC. As 
can be seen, there is a remarkable effect of numerical resolution on the 
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Cummins Engine Data 
Cylinder Bore 
Stroke 

Compression ratio 
Displacement 

Number of spray nozzle orifices 
Nozzle hole diameter 
Spray axis, angle from head 
Combustion chamber 
Piston crown 
Engine speed 
Overall equivalence ratio 
Air flow rate 
Fuel flow rate 
Intake tank pressure 
Intake tank temperature 


139.7 mm 
152.4 mm 
13.23 

2.33 liters 
8 

0.2 mm 
18° 

Quiescent 
Mexican hat type 
1500 r/min (constant) 
0.6 

0.353e-2 kg/cycle 
0.144e-3 kg/cycle 
148.2 kPa 
302 K 


Engine temperatures at IVC (constant during the cycle) 


Cylinder wall 
Cylinder head 
Piston surface 
Valves face 

Mass average gas temperature at IVC 
Cylinder pressure at IVC 
Swirl number 
Fuel 

Injection starts 
Injection ends 


405 K 
486 K 
578 K 
773 K 
359 K 
157.9 kPa 
1.0 


Tetradecane 
18° btdc 
11° atdc 




predicted spray penetration at TDC. It should be noted that these 
computations were made using the atomization model of Reitz [18] with no 
modification to the vaporization routines. 

Figure 8 shows computations that include the effect of combustion, made 
using the same grids as those for the results in Figs. 6 and 7. As can be seen, 
the fine grid computation gives a higher predicted peak cylinder pressure 
than the coarse grid computation. However, the predicted pressures still 
underpredict the measured pressures substantially. The results shown in Fig. 
9 demonstrate that the details of the combustion model have a large effect on 
the predicted cylinder pressures, particularly for the fine grid computations. 
For Fig. 9 the computations used a reduced value of the pre-exponential 
constant, viz. Ki = 7.68xl0 8 . This value of Ki was found to give the highest 
possible peak cylinder pressure for both the fine and the coarse grid cases. 
Further reductions in Ki caused unrealistically long ignition delays. It is seen 
that the discrepancy between measurements and predicted cylinder pressures 
is reduced significantly for the fine grid computation with the revised 
combustion model constant. 



Figure 6 Plan view of sprays showing effect of grid resolution on predicted 
spray penetration. TDC, 1500 rev/min, computations made without 
considering combustion. Top: Coarse grid which uses 25x6x18 cells. Bottom: 
Fine grid which uses 24x12x18 cells with non-uniform azimuthal spacing. 
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Figure 7 Elevation view showing effect of grid resolution on predicted spray 
penetration for the conditions of Fig. 17. Top: Coarse grid which uses 25x6x18 
cells. Bottom: Fine grid which uses 24x12x18 cells with non-uniform 
azimuthal spacing. 
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The reason for the improved pressure predictions is seen in Fig. 10. Figure 10 
shows the predicted temperature contours and spray drop locations in the 
chamber at TDC in a plane through the center of the spray for the coarse and 
fine grid cases of Figs. 6, 7 and 9 (i.e., using the revised combustion model 
constant Kj). As can be seen in the fine grid results, the high temperature 
contours are now located closer to the edge of the bowl, in better agreement 
with the experimental observations of Yan and Borman [21] who, using a 
radiation probe mounted in the cylinder head, found that the flame reached 
the piston bowl outer surface between 6 deg BTDC and 1 deg ATDC. However, 
the coarse grid temperature contours show that the flame is still confined to 
the region near the nozzle. (The position of the high temperature regions in 
the combustion cases of Fig. 10 is seen to be correlated with the spray 
penetrations (without combustion) shown in Fig. 7.) 

The sensitivity of the fine grid results to the combustion model constants 
indicates that further optimization of the results in Fig. 9 will require 
improved combustion models. For example, the simplified Arrhenius-type 
combustion model used in this study does not specifically account for ignition 
processes nor does it include the effect of turbulent mixing that is known to 
be important in diesel combustion. This is the subject of Task 5 below. 

The atomization model itself has been improved by extending the stability 
analysis to consider an evaporating liquid surface. In this analysis [24] the 
evaporation flux of a perturbed surface with local curvature R was assumed 
to be proportional to that of a spherical droplet of radius R, and the kinematic 
condition (jump condition) at the liquid-gas interface is modified to include 
the effect of surface evaporation velocity. The recoil force due to evaporation 
is also taken into account in the interface normal stress equation [25,26]. 

A new dispersion equation was obtained as shown in Fig. 11, where Va is a 
dimensionless parameter representing the evaporation level. Figure 11 
shows that surface evaporation reduces both the growth rate and 
wavenumber of the most unstable waves; hence the surfaces waves break 
slower and produce larger droplets. The influence of this phenomenon on 
the atomization process is being evaluated. 

An additional consideration is the influence of the compressibility of the 
ambient gas on the atomization process. In the stability analysis of Reitz [18], 
it was assumed that both liquid and gas are incompressible. The validity of 
this assumption is questionable for the high injection pressures of interest in 
diesel engines since the Mach number often exceeds 0.3. To study the gas 
compressibility on high speed spray, we have adopted a method where the 
dynamic pressure in the gas includes the effect of Mach number following 
Bradley [27,28]. The effect of gas compressibility on atomization is being 
assessed. 


10000 



Crank Angle (degrees) 

Figure 8 Predicted cylinder pressure using the unmodified combustion 
model with Ki = 3.68x1 0 10 showing effect of grid resolution for the coarse and 
fine grids of 
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Figure 9 Predicted cylinder pressure using the modified combustion model 
with Ki = 7.68x10 s showing effect of grid resolution for the coarse and fine 
grids of Fig. 6. 






Figure 10 Elevation view showing effect of grid resolution on predicted spray 
penetration and temperature contours for the conditions of Fig. 6. Top: 
Coarse grid which uses 25x6x18 cells (h=2660 K, 1=917 K). Bottom: Fine grid 
which uses 24x12x18 cells with non-uniform azimuthal spacing (h=2720 K, 
1=1090 K). 
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Figure 11 Predicted wave growth rate versus wavenumber for different 
evaporation levels from the modified dispersion relationship. 

3. HC and NOx emission 

The combustion model specifies the rate of conversion from reactants to 
products using a characteristic time which is controlled by the longer of the local 
turbulent mixing or laminar kinetics times and is formulated such that correct 
thermodynamic equilibrium is recovered in the burned gas [17]. The correct 
burned-gas temperature is needed for accurate NOx predictions using the 
Zeldo'vich kinetics model which has been included in KIVA . The combustion 
and NOx models are included in Subroutine Chmprn in Appendix 3. 

The combustion model automatically predicts unbumed hydrocarbon 
(HC) in regions of low gas temperature or high equivalence ratio where the 
characteristic chemistry time is large and the conversion of fuel to products is 
arrested. Unburned HC and NOx predictions will be examined in the model 
application project (see Table 2, page 33). 

4. Multicomponent Fuel Vaporization 

Our initial attempts to model diesel engine combustion phenomena revealed 
major qualitative discrepancies between computations and experiments [19] 

(see also Task 2 above). These discrepancies motivated a critical evaluation of 
single component drop vaporization models. 

The fuel considered in the study was hexadecane (Tcr=725 K), and a single 50- 
pm diameter droplet was injected at a velocity of 100 m/ s into a cylinder 
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containing air at 40 atm and 800 K. The mesh used in this phase had 20x1x60 
cells (1x2 mm cells). 

The original KIVA vaporization model is a low pressure model that is based 
on the following assumptions : a. spherical symmetry, b. single component 
fuel, c. uniform drop temperature, d. constant fuel density and, e. unity Lewis 
number. Assumptions d. and e. have been removed in the present effort, and 
the details are summarized below. The preparatory phases of modeling 
multi-component fuels (assumptions b. and c.) has been completed. The 
spherical symmetry assumption was considered to be reasonable for stable, 
non-breaking drops. 

Since the injected drops undergo a significant temperature variation during 
their life-time and the properties of the fuel change considerably, the 
governing equations for vaporizing drops were rewritten to include variable 
liquid density which was not accounted for in the original KIVA. 

Figure 12 shows the predicted drop size variation (drop radius non- 
dimensionalized by initial radius) as a function of time using the original 
KTVA vaporization model, and using the new model that accounts for 
variable liquid density. With variable density, the drop size increases in the 
heat-up period. This size increase is accompanied by an increase in the 
available surface area which increases the vaporization rate resulting in a 
shorter drop lifetime. 

An additional factor that was considered is the influence of mass transfer on 
the heat transfer for vaporizing drops. Various corrections have been 
proposed in the literature but the one used in KTVA assumes that the Lewis 
number of the fuel vapor is unity. Computational results under the 
conditions of interest show that the Lewis number is typically between 4 and 
5. This observation has motivated the use of the correction of Priem et al. [29] 
in place of the KIVA correction. 

As can be seen in Fig. 12, the variable Lewis number effect (which represents 
superheating of the vapor around the drop) also leads to shorter drop 
lifetimes. The modified version of KTVA now includes both variable fuel 
density and Lewis number effects, and the results are seen in Fig. 12 to lead to 
a significant reduction in the predicted drop lifetime. Other results show that 
these reductions in lifetime are not as large at lower values of the drop-gas 
relative velocity, as would be expected from the form of the correction terms. 

Work is also in progress on implementing a multi-component fuel 
vaporization model into KIVA that is based on the works of Refs. [30,31]. 
Preliminary results are shown in Fig. 13 for a 50 jim diameter, 50% 
pentane/50% octane drop at 300 K injected into a 800K, 4 MPa nitrogen 
environment. The drop was injected at 100 m/s. 
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r/r for various models 



Figure 12 Predicted drop size variation with time showing the effect of 
variable fuel density (var p) and non-unity Lewis number (z only). Modified 
KIVA results includes both density and Lewis number effects. 



Figure 13 Predicted fuel drop location and concentration contours at 2 ms 
after injection (the injector is located at the top center of the diagrams). 50 nm 
diameter, 50% pentane/50% octane drop at 300 K injected into a 800K, 4 MPa 
nitrogen environment. Chamber is 40 mm in diameter. 
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5. Diesel ignition and combustion 

Our preliminary modeling of diesel ignition and combustion used kinetic 
constants from the ignition study of Bergeron and Hallet [32] to specify the 
Arrhenius rates for the laminar chemistry time in the combustion model 
[19,23]. In the current effort, the laminar and turbulent characteristic time 
combustion model of Abraham et al. [33] that was used for the spark-ignited 
engine studies (see Task 2, above) has been extended to allow predictions of 
ignition [34]. 

The homogeneous charge, compression-ignition engine experiments of Boggs 
[35] were used for this study. These experiments were performed in the CFR 
engine at the Center using ethylene fuel. Six cases were studied that include 
two engine speeds (600 and 1200 rev/ min), two swirl ratios (1.2 and 8) and two 
compression ratios (10.5 and 6.3 - equivalence ratios of 0.24 and 0.4). Details of 
the engine geometry are given below. The residual amounts and exhaust 
temperature were estimated by using the method of Yuh and Mirsky [36]. 

CFR Engine specifications 


cycle 

Bore 

Stroke 

Connecting Rod Length 
Displacement 
Valve 
Lift 

Intake timing 
Exhaust timing 
Cooling System 


4-Stroke 
83.1 mm 
114.3 mm 
254.0 mm 
620 c.c. 

Supercharge cam 
7.93 mm 

Open 15 0 BTC, Close 50 0 ABC 
Open 50 0 BBC, Close 15 0 ATC 
Evaporative, atmospheric pressure 
Water coolant 


The time rate of change of the partial density of species i , due to conversion 
from one chemical species to another, is given by dYj/dt = -( Yi - Yf)/tc [33], 
where Yi is the mass fraction of species i , and the * indicates local and 
instantaneous thermodynamic equilibrium values. t c is the characteristic 
time for reaction which is assumed to be the sum of the laminar (high 
temperature) chemistry time , the turbulence mixing time tf , and ignition 
(low temperature) chemistry time ti, which was not present in the original 
combustion model, i.e., tc = + tt + ti . 

The ignition characteristic time ti was modeled using data from elementary 
initiation reactions of ethylene [37] and has the Arrhenius form. Since this 
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timescale represents many possible reaction paths, it was felt justifiable to 
adjust the pre-exponential constant slightly to give the best overall agreement 
with the experiments. It was found to be possible to match all cases 
reasonably well with one set of ignition model constants. Crevice flow effects 
were found to be significant for this engine, and the crevice flow model of 
Reitz and Kuo [38] was implemented in KTVA. It was necessary to measure 
the piston-cylinder-liner crevice volumes on the engine. 

Typical results showing comparisons between measured and predicted 
cylinder pressures are given in Fig. 14. The level of agreement is quite good. 

It should be noted that this level of agreement was not possible if the 
turbulent mixing timescale was not included as a parameter in the 
combustion model. 


0608 p.f.plot 



Figure 14 Comparison between computed and measured cylinder pressures 
for compression ignited, homogeneous charge engine combustion (engine 
speed 1200 rev/min, compression ratio 11.5, swirl ratio 8.0, equivalence ratio 
0.24). 

This important conclusion implies that thermal mixing phenomena control 
the rate of combustion, even in this engine that is characterized by the 
absence of a propagating flame. The characteristic timescales are shown in 
Fig. 15. As can be seen, ignition is controlled at low temperatures by the 
ignition timescale. Notice that the high temperature laminar chemistry 
timescale never plays a role in determining the combustion rate since it is 
always smaller than the other characteristic times. 
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Figure 15 Characteristic combustion times versus cylinder gas temperature 
for the compression-ignited, homogeneous-charge engine case of Fig. 14. 

The location of the combustion regime for this engine has been compared 
with the internal combustion engine regimes proposed by Abraham et al. [39] 
which consider the Damkohler number (ratio of turbulence mixing time to 
the reaction time) and the turbulent Reynolds number (based on the integral 
length scale). The present engine combustion was found to be located in the 
distributed reaction regime once combustion starts. 

The corresponding wall heat flux predictions are compared with the 
experimentally measured values in Fig. 16. The phasing of the heat flux is 
well predicted, although the peak values are under-estimated by amounts 
that are similar to those for the spark-ignited engine shown in Fig. 1. The 
implication of these results about the necessity to include the effect of 
combustion in wall heat transfer models is being investigated. 

The next phase of the work will involve the above characteristic time 
ignition model, and other ignition models to improve spray combustion 
predictions (see Task 2 above). 

6. Soot and Radiation 

Work will begin on this effort in the next year (see Future Work). 


26 




0608 q.f.plot 



Figure 16 Comparison between predicted and measured wall heat flux for 
compression ignited, homogeneous charge engine combustion (engine speed 
1200 rev/min, compression ratio 11.5, swirl ratio 8.0, equivalence ratio 0.24). 

7. Intake Flow Modeling 

The intake flow work is proceeding on schedule. Initial phase one 
calculations are successfully being performed to evaluate simple flow fields 
and different grid schemes. In addition, grid generation for realistic multi- 
valve port geometries is nearing completion. 

Intake flow modeling is being implemented using KIVA-3. This is a new, 
unreleased version of KIVA for which the Engine Research Center is serving 
as a beta test site. The major advancement of KIVA-3 over previous versions 
is the use of unstructured grids. This provides an efficient means of 
representing complex, multi-domain geometries such as the intake port, 
valve, and cylinder. In addition, KIVA-3 retains the ability to model moving 
grids and allows different regions to become disconnected during a 
calculation. 

KIVA-3 is now tested and running on the Cray machines used at the ERC. 

An interface between KIVA-3 and the graphical post-processing program, 
PLOT3D, has also been implemented. Initial calculations on simple (but 
multi-domain) geometries are being carried out to test feasibility of gridding 
schemes and for initial examination of intake flow fields. 

An example calculation is shown in Fig. 17. The initial gridding and sample 
calculation used a stationary valve head and stem with a moving piston and 
considered the intake stroke. The 3-D view of the velocity field in Fig. 17 was 
obtained by converting KTVA-3's post-processor to give PLOT3D compatible 
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Figure 17 Intake flow computation using KIVA-3 for a centrally located 
stationary valve at 120 degrees after top-dead-center. 
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output. The computational grid was created using the KTVA-3's pre- 
processor. A constant pressure boundary condition was used at the top of the 
intake port. The flow enters the cylinder and is deflected radially outward by 
the valve head. The cylinder wall deflects the flow downwards and a large 
scale recirculation pattern is established. Despite the coarse grid used in these 
initial computations, the flow patterns are very similar to visualization 
experiments by Eaton and Reynolds [40] in a motored engine with a centrally 
located valve. 

Grid generation remains a major task for the intake flow calculations, 
since it is so time consuming to create a complex 3D grid. Work is 
underway to replace KIVA-3's simple grid generator with the GRIDGEN 
code developed by General Dynamics under Air Force contracts. This code 
has been selected over the EAGLE code because it is more flexible and it 
makes use of our Silicon Graphics IRIS workstation. Also, GRIDGEN is 
being used by many other users (we have joined a grid generation users 
group (NAS SigGrid, NASA Ames)). 

Grid generation with GRIDGEN is a three step process. The first step is the 
division of the geometry into logical blocks. The second step is the formation 
of a computational grid on each face of each block. This establishes the 
surface grid that will be used to generate the volume grid in the third step. 

The final step is the generation of the 3D volume grid. This is a fairly 
intensive calculation that is currently implemented only on CRAY 
computers running COS. Modifications to run under TJNICOS or on the 
Stardent Titan are currently being investigated. Results of this step are 
transferred back to the Iris for visualization. The file containing the 3D grid 
information is then reformatted for input to KIVA-3. 

The basic geometry used in step one of this process is obtained from digitized 
data of the intake manifold. An expanding flexible foam was used to 
transform the manifold's internal features into a mold form. Liquid acrylic 
was poured around the foam form to create a full size replica of the original 
manifold geometry. The acrylic casting was then cut into slices. Features 
machined into the acrylic block (see Fig. 18) allow the longitudinal position of 
each slice to be determined. 

A computerized video measurement technique transposes the physical 
geometry contained in the cross sections into coordinate data. This technique 
digitizes the images using edge finding algorithms, and calculates the line 
coordinates using parametric splines. The output is the three dimensional 
coordinates that correspond to the cross section geometry. The technique is 
repeated for each slice until all geometry is recorded. Finally, the cross section 
geometry is combined and formatted to serve as input for step 1 of the grid 
generation process. 
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Figure 18 Schematic showing foam mold of intake port and acrylic block 
used to obtain digitized geometry. 

Figure 19 shows a grid of a Caterpillar engine intake port (the interior grid is 
removed for clarity) that has been generated using this method. 


Various approaches for implementing a moving valve into KIVA-3 are 
currently being evaluated. The moving valve requires that different regions 
of the grid move relative to each other. This will implemented with a 
variation of the "snapper" algorithm used to move the piston in KIVA-3. 

This algorithm maintains grid point connectivity by snapping grid points 
between neighboring cells as the grid moves. The snap occurs when the gird 
has moved a specified percentage of a cell distance. A method of using the 
snapper algorithm with the moving valve has been developed and is shown 
schematically in Fig. 20. Figures 20 (a) and (b) show the grid as valve begins to 
move upwards; all vertices are active in the calculation until the grid reaches 
a deformation limit. At the deformation limit (Fig. 20 (c) ), the snapper 
algorithm inactivates the six open circle vertices and moves the eight square 
vertices down to the valve surface. The six vertices on the valve face are 
moved down into the cylinder region and replaced by new vertices. Volume 
weighted interpolation is used to reassign all properties associated with 
snapped vertices and their cells. 
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Figure 19 Grid generation for a single cylinder Caterpillar engine intake port. 
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Figure 20 Schematic showing use of "snapper" algorithm for moving valve. 
Open circles and squares are used to show the movement of vertices as the 
valve closes. 
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B. Model Application 

The model application tasks in Phase B are presented in Table 2. The 
individual projects are discussed below in the Future Work section. 


Table 2 - Model Application 


Year 

Topic / Study 

9/1/89 

9/1/90 

1 

9/1/90 

9/1/91 

2 

9/1/91 

9/1/92 

3 

9/1/92 

9/1/93 

4 

9/1/93 

9/1/94 

5 

1 

Ignition delay 


X 

X 

X 


2 

Fuel-air mixing 



X 

X 

X 

3 

Wall impingement 



X 

X 

X 


FUTURE WORK 

The proposed work schedule for the period 6/91 - 9/92 for Parts A (Model 
Implementation) and B (Model Application) can be seen in Tables 1 and 2. 

The plan of work for Part C (Engine Experiments) is discussed on page 35. 
Details of some of the individual tasks are described below. 

Part A - Multi-component fuel vaporization and atomization 

We will continue to study the influence of vaporization on atomization [24] 
and assess the effect on diesel combustion modeling. In addition we will 
extend our implementation of the multi-component fuel vaporization model 
in KIVA. To save computer time and storage, we are exploring methods 
where the droplet internal temperature and concentration fields can be 
solved using at most ten internal grid points per drop. As the drop vaporizes, 
the number of internal grid points is reduced and eventually, when the drop 
size becomes sufficiently small, the lumped parameter approach of KIVA is 
recovered. 

Part A - Diesel Ignition 

We will extend our ignition model to consider spray ignition. In this case the 
high temperature constant volume bomb experiments of Kamimoto [20] will 
be used as a starting point, since they also include cases with combustion. Our 
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previous diesel combustion computations [19,23] that modeled the data of 
Yan and Borman [21] for a Cummins NH engine will be continued. We will 
also be making use of measured pressure data from the Caterpillar engine 
experiments (see Part C). As part of the ignition modeling effort, we will 
explore the use of the Theobald ignition model [41]. 

Part A - Soot and Radiation 

Soot model development will start by considering the soot model proposed by 
Gentry et al. [42] for coal slurry applications which will be implemented. Due 
to the complexity of the soot formation and destruction processes, it is 
anticipated that soot concentration predictions will be qualitative only. It is 
planned that the radiation model of Chang and Rhee [43] will be adapted for 
use in this study. This model assumes the Rayleigh-limit for the soot 
emissivity and a band model for the gas emissivity. More complicated 
radiation models that require solution of the radiation transport equation and 
that account for scattering by the spray drops are available [44]. However, in 
view of the uncertainties about soot concentrations, this additional 
complexity is not considered to be warranted at the present time. Much work 
is being done in this area and it will be monitored for its appropriateness for 
inclusion in KIVA. 

Part A - Intake flow modeling 

We expect to be able to perform computations for realistic intake port and 
valve geometries in the near future. Comparisons will be made to steady 
state flow bench experimental data. Swirl measurements and flow 
visualizations will be made for comparison to the KIVA calculations using 
the Caterpillar multi-valve intake port. Any discrepancies between 
experiments and calculations will be investigated. This may require 
improvements to the turbulence model in KIVA so that the high shear 
through the valve curtain is more accurately modeled. 

Part B - Ignition delay 

Computations of the influence of diesel ignition delay on emissions will be 
made. The data of Yan and Borman [21] will be used to verify the accuracy of 
the model. The goal is to explore research issues such as the influence of the 
injection timing and the number of injector holes on ignition delay. Over- 
mixed vaporized fuel is thought to be responsible for the increase in HC 
emissions found with an increased ignition delay. The influence of the fuel 
injection parameters on multi-component fuel vaporization and mixing 
during the ignition delay and the early stages of the diffusion burning phase 
will be studied in later phases of the project. 
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Part B - Fuel-air mixing, Spray/wall impingement 

Computations of fuel-air mixing and spray wall impingement will begin in 
the summer of 1992. The location and extent of over- and under-mixed fuel- 
air regions within the combustion chamber influence engine operation and 
omissions. Particulates are thought to be formed when under-mixed regions 
are exposed to high temperatures. The influence of the fuel injection 
parameters and in-cylinder flow details on mixing and combustion during 
the diffusion burning phase will be assessed. 

The effect of wall impingement on in-cylinder fuel-air distributions and 
combustion will also be considered. In particular, differences between sprays 
whose liquid core extends to the wall and sprays whose atomization is 
completed before wall impingement will be assessed. This is of interest 
because of the current trends toward smaller engines and higher injection 
pressures, both of which promote increased wall impingement. Wall 
impingement processes are thought to influence NOx and soot emissions and 
wall heat transfer. This will be examined in the study. 

PART C 

An experimental program has been initiated at the University of Wisconsin 
with funding from Caterpillar and from DOE/ NASA-Lewis for engine 
experiments to generate data needed for KIVA validation. The work will 
comprise three tasks: 1. in-cylinder gas velocity and turbulence 
measurements, 2. combustion visualization experiments and 3. 
measurements of engine-out emissions. 

In setting up the experimental program, effort has been taken to ensure that 
the modifications to the engine needed for optical access will be minimal so 
that the results will be representative of actual engines. The experiments will 
be conducted using dual and single intake valve configurations and with 
standard and re-entrant piston designs. 

The in-cylinder velocity measurements will be made using Particle Image 
Velocimetry [45,46]. The method has already been used at the University of 
Wisconsin [47], and techniques for rapid data analysis are available at the ERC 
[48]. 

Combustion visualization and photography experiments will be conducted 
through a window that replaces one of the exhaust valves. If window fouling 
problems prove to be insurmountable, a radiation probe will be used that has 
been developed and used successfully at the ERC [49]. Engine-out particulate 
measurements will be made using a mini-dilution tunnel sized for the CAT 
engine. It is anticipated that these emission measurements will be useful in 
assessing the validity of KIVA’s soot models. 


35 



ACKNOWLEDGEMENTS 


Progress on the grant has benefited from numerous interactions with faculty 
and students working on related research projects at the Engine Research 
Center under ARO and other funding. The donation of computer time by 
CRAY at Mendota heights and at the San Diego Supercomputer Centers are 
also appreciated. 

REFERENCES 

1. Amsden, A. A., Ramshaw, J. D., O’Rourke, P. J. and Dukowicz, J. K., 

"KIVA: A Computer Program for Two- and Three-Dimensional Fluid Flows 
with Chemical Reactions and Fuel Sprays," Los Alamos Scientific 
Laboratory Report LA-1 0245-MS, February 1985. 

2. Amsden, A. A., Butler, T. D. and O'Rourke, P. J. "The KTVA-II Computer 
Program for Transient Multi-Dimensional Chemically Reactive Flows with 
Sprays," Society of Automotive Engineers Technical Paper 872072, 1987. 

3. Glikin, P.E. "Fuel Injection in Diesel Engines," Proceedings of the 
Institution of Mechanical Engineers, Vol. 199, No. 78, pp. 1-14, 1985. 

7. Greeves, G. "Response of Diesel Combustion to Increase in Fuel Injection 
Rate," Society of Automotive Engineers Technical Paper 790037, 1979. 

8. Greeves, G. and Wang, C.H.T. "Origins of Diesel Particulate Mass 
Emissions," Society of Automotive Engineers Technical Paper 810260, 1981. 

6. Yoshikawa, S., Furusawa, R., Arai, M. and Hiroyasu, H. "Optimizing spray 
behavior to improve engine performance and to Reduce Exhaust Emissions 
in a Small DI Diesel Engine" Society of Automotive Engineers Technical 
Paper 89046, 1989. 

7. Cartellieri, W.P. and Herzog, P.L. "Swirl supported or Quiescent 
Combustion for 1990s Heavy-Duty DI Diesel Engines - An Analysis," Society 
of Automotive Engineers Technical Paper 880342, 1988. 

8. Siegla, D.C. and Alkidas, A.C. "Evaluation of the Potential of a Low-Heat- 
Rejection Diesel Engine to Meet Future EPA Heavy-Duty Emission 
Standards," Society of Automotive Engineers Technical Paper 890291, 1989. 

9. Gonzalez D., M.A., Borman, G.L. and Reitz, R.D. "A Study of Diesel Cold 
Starting using both Cycle Analysis and Multidimensional Calculations," 
Society of Automotive Engineers Technical Paper 910180, 1991. 

10. Abraham, J. and Bracco, F.V. "Fuel-Air Mixing and Distribution in a 
Direct-Injection Stratified-Charge Rotary Engine," Society of Automotive 
Engineers Technical Paper 890329, 1989. 

11. Kuo, T.-W. and Reitz, R.D. "Computations of Premixed-Charge 
Combustion in Pancake and Pent-roof Engines," Society of Automotive 
Engineers Technical Paper 890670, 1989. 

12. Reitz, R.D. "Three-Dimensional Modeling of Diesel Spray Combustion 
and Emissions," Proceedings of the DOE Automotive Technology 


36 



Development Contractors Coordination Meeting, Dearborn, MI, October 22- 
25, 1990, SAE P-243. 

13. Reitz, R.D., and Rutland, C.J. "Three-Dimensional Modeling of Diesel 
Engine Intake Flow, Combustion and Emissions," Proceedings of the DOE 
Automotive Technology Development Contractors Coordination Meeting, 
Dearborn, MI, October 28-31, 1991. 

14. Reitz, R.D., and Rutland, C.J. "Three-Dimensional Modeling of Diesel 
Engine Intake Flow, Combustion and Emissions," Society of Automotive 
Engineers Paper 911789, 1991. 

15. Naber, J.D. and Reitz, R.D. "Modeling Engine Spray /Wall Impingement" 
Society of Automotive Engineers Paper 880107, 1988. 

16. Gonzalez D., M.A., Borman, G.L. and Reitz, R.D. "A Study of Diesel Cold 
Starting using both Cycle Analysis and Multidimensional Calculations," 
SAE paper 910180, 1991. 

17. Reitz, R.D., 'Assessment of Wall Heat Transfer Models for Premixed- 

Charge Engine Combustion Computations,’ Submitted to Society of 

Automotive Engineers Transactions, SAE paper 910267, 1991. 

18. Reitz, R.D. "Modeling Atomization Processes in High-Pressure Vaporizing 
Sprays," Atomisation and Spray Technology, 3, pp. 309-337, 1987. 

19. Gonzalez D., M.A. and Reitz, R.D. "Modeling Diesel Engine Spray 
Vaporization and Combustion," Presented at the First International KIVA 
Users Group Meeting, Detroit, MI, February 24, 1991, and at the Fifth 
International Conference on Liquid Atomization and Spray Systems, NIST 
Gaithersburg, MD, USA, July 15-18, 1991. 

20. Kamimoto, T., Yokota, H. and Kobayashi, H. "Effect of High Pressure 
Injection on Soot Formation Processes in a Rapid Compression Machine to 
Simulate Diesel Flames," Society of Automotive Engineers Paper 871610, 
1987. 

21. Yan, J. and Borman, G.L. "Analysis and in-cylinder measurement 
of particulate radiant emissions and temperature in a direct 
injection diesel engine," SAE Paper 881315 (1988). 

22. Lei, N. "A cycle simulation program for the dynamic operation of 
a single cylinder direct injection diesel," MS Thesis, Department of 
Mechanical Engineering, University of Wisconsin-Madison (1988). 

23. M. A. Gonzalez D., Z. W. Lian and R. D. Reitz, "Modeling Diesel Engine 
Spray Vaporization and Combustion," Submitted for publication for the 1992 
SAE Congress and Exposition, Detroit, MI. 

24. Z. W. Lian and R. D. Reitz, "Effect of Vaporization and Gas 
Compressibility on Liquid Jet Breakup," Submitted for Publication, Physics 
of Huids, 1991. 

25. Prosperetti, A. and Plesset, M. S., 'The Stability of An Evaporating Liquid 
Surface,' Phys. Fluids 27, 7, 1984. pp. 1590-1602. 

26. Hiquera, F. J., 'The Hydrodynamic Stability of An Evaporating Liquid,' 
Phys. Fluids 30, 3, 1987. pp. 679-686. 

27. Bradley, D. 'On the Atomization of Liquid By High-Velocity Gases,' 

J. Phys. D: Appl. Phys., Vol. 6, 1973. 


37 


28. Bradley, D. 'On the Atomization of A Liquid by High-Velocity Gases: II,' 

J. Phys. D: Appl. Phys., Vol. 8, 1975. 

29. Priem, R. J., Borman, G. L., El-Wakil, M. M., Uyehara, O. A., and Myers, P. 
S., "Experimental and Calculated Histories of Vaporizing Fuel Drops", 
NACA TN 3988, 1957. 

30. Jin, J. D., and Borman, G. L., "A Model for Multicomponent Droplet 
Vaporization at High Ambient Pressures", Society of Automotive Engineers 
paper 850264, 1985. 

31. Abramazon, B., and Sirignano, W. A., "Approximate Theory of a Single 
Droplet Vaporization in a Convective Field", ASME/JSME Thermal 
Engineering Conference, Vol. 1, 1987, pp. 11-18. 

32. Bergeron, C.A. and Hallett, W.L.H. "Autoignition of Single Droplets of 
Two-Component Liquid Fuels," Combust. Sci. and Tech., 65, pp. 181-194, 
1989. 

33. AbrahamJ., Bracco,F.V., and Reitz,R.D., "Comparison of Computed and 
Measured Premixed Charged Engine Combustion", Combustion and Flame, 
60, 1985. 

34. Kong, S.-C., Ayoub, N., and Reitz, R.D., "Modeling Combustion in 
Compression Ignition Homogeneous Charge Engines," Submitted for 
publication for the 1992 SAE Congress and Exposition, Detroit, MI. 

35. Boggs,D.E., Ph.D. Thesis, University of Wisconsin-Madison,1990. 

36. Yuh,H.J., and Mirsky,W., "Schlieren-Streak Measurements of 
Instantaneous Exhaust Gas Velocity from a Spark Ignition Engine", Society 
of Automotive Engineers Paper 741015, 1974 

37. Westbrook,C.K., and Dryer,F.L., Combust. Sci. Technol., 20, 125, 1979. 

38. Reitz, R.D. and Kuo, T.-W. "Modeling of HC Emissions Due to Crevice 
Flows in Premixed-Charge Engines," Society of Automotive Engineers 
Paper 892085, 1989. 

39. AbrahamJ., William, F.A., and Bracco,F.V.,"A discussion of turbulent 
flame structure in premixed charges", Society of Automotive Engineers 
Paper 850345, 1985 

40. Eaton, A.R. & Reynolds, W.C., "Flow structure and mixing in a motored 
axisymmetric engine," Society of Automotive Engineers Paper 890321, 1989. 

41. Theobald, M.A. and Cheng, W.K. "A Numerical Study of Diesel Ignition," 
American Society of Mechanical Engineers Paper 87-FE-2, 1987. 

42. Gentry, R.A., Daly, B.J. and Amsden, A.A. "KIVA-COAL: A Modified 
Version of the KTVA Program for Calculating the Combustion Dynamics of 
a Coal-Water Slurry in a Diesel Engine Cylinder," Los Alamos National 
Laboratory Report LA-1 1045-MS, August 1987. 

43. Chang, S.L. and Rhee, K.T. "Computation of Radiation Heat Transfer in 
Diesel Combustion," Society of Automotive Engineers Paper 831332, 1983. 

44. Menguc, M.P., Viskanta, R. and Ferguson, C.R. "Multidimensional 
modeling of Radiative Heat Transfer in a Diesel Engine," Society of 
Automotive Engineers Paper 850503, 1985. 

45. Reuss, D.L., Adrian, R.J., Landreth, C.C., French, D.T. and Fansler, T.D. 
"Instantaneous Planar Measurements of Velocity and Large Scale Vorticity 


38 


and Strain Rate in an Engine Using Particle-Image Velocimetry," SAE 
Technical Paper 890616, 1989. 

46. Reuss, D.L., Bardsley, M., Felton, P.G., Landreth, C.C. and Adrian, R.J. 
"Velocity, Vorticity, and Strain-Rate Ahead of a Flame Measured in an 
Engine Using Particle Image Velocimetry," SAE Technical Paper 900053, 
1990. 

47. Ghandi, J.B. 'Velocity Field Characteristics in a Motored Two-Stroke 
Ported Engine,' PhD Thesis, University of Wisconsin-Madison, 1991. See 
also ERC Newsletter, Vol. 4, No. 1, April 1991. 

48. Goetsch, D. and Farrell, P.V. 'Optical LCTV Correlator for Particle Image 
Velocimetry,' Optics Letters, 14, 978, (1989). 

49. Yan, J. "Analysis and In-cylinder Measurements of Local and 
Hemispherical Particulate Radiant Emissions and Temperatures in a Direct 
Injection Diesel Engine," Ph.D. Thesis, University of Wisconsin, 1988. 


39 


APPENDIX 1 KIVA USERS GROUP ACTIVITY 


The KIVA computer code and its predecessor codes were written at the Los 
Alamos National Laboratories, New Mexico [1,2]. These codes have seen 
increasing use by university researchers, government laboratories and engine 
industries in the U.S. and abroad since the early 1970s. A U.S. KIVA Users 
group was formed during the DOE Diesel Working Group Meeting at the 
University of Wisconsin-Madison in the Fall of 1989 . The purpose of this 
Users Group is to promote the use of KTVA and to facilitate exchange of 
information among KIVA users. Since that time KTVA Users Groups have 
also been formed in Europe and Japan. 

The U.S. KTVA Users Group currently numbers about 80 organizations. We 
have served as the editors of the KIVA Users Newsletter and four issues have 
already been distributed to the membership (copies of the newsletters are 
attached in this Appendix). We also co-organized the First International 
KTVA Users Group Meeting with CRAY Research which took place in 
February 1991 in Detroit, MI. This meeting was attended by 75 registrants. A 
similar meeting is currently being planned for 1992. Details will be released 
in the 5th KIVA Users Newsletter later this month. 
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APPENDIX 2 LIST OF RELATED PUBLICATIONS AND ABSTRACTS 
(papers 1 and 4-8 partially funded under grant NAG 3-1087) 

1. Reitz, R.D. "Assessment of Wall Heat Transfer Models for Premixed-Charge 
Engine Combustion Computations," Accepted for publication in SAE 
Transactions, SAE Paper 910267. 

Two-dimensional computations of premixed-charge engine combustion were made using the 
KTVA-II code. The purpose of the study was to assess the influence of heat transfer and turbulence 
model boundary conditions on engine combustion predictions. Combustion was modeled using a 
laminar- and turbulent-characteristic-time model. Flow through the piston-cylinder-ring crevice 
was accounted for using a phenomenological crevice-flow model. The predictions were compared 
to existing cylinder pressure and wall heat transfer experimental data under motoring and fired 
conditions, at two engine speeds. Two different wall heat transfer model formulations were 
considered. The first is the standard wall function method. The second is based on solutions to 
the one-dimensional unsteady energy equation, formulated such that the standard wall function 
method is recovered in the quasi-steady limit. Turbulence was modeled using the standard k-e 
turbulence model equations. However, the turbulence model boundary conditions were modified to 
account for compressibility effects by using a coordinate transformation in the wall region. The 
results show that the details of wall heat transfer and turbulence model boundary conditions 
influence heat transfer predictions greatly through their influence on the flame speed and the 
flame structure in the vicinity of the wall. Inclusion of compressibility and unsteadiness effects 
leads to increased wall heat flux values that agree better with measurements. 

2. Epstein, P., Reitz, R.D. and Foster, D. "Computations of Two-Stroke Engine 
Cylinder and Port Scavenging Flows," Accepted for publication in SAE 
Transactions, SAE paper 910672. 

A modification of the computational fluid dynamics code KIVA-II is presented that allows 
computations to be made in complex engine geometries. An example application is given in 
which three versions of KTVA-II are run simultaneously. Each version considers a separate block 
of the computational domain, and the blocks exchange boundary condition information with 
each other at their common interfaces. The use of separate blocks permits the connectedness of 
the overall computational domain to change with time. The scavenging flow in the cylinder, 
transfer pipes (ports), and exhaust pipe of a ported two-stroke engine with a moving piston was 
modeled in this way. Results are presented for three engine designs that differ only in the angle 
of their boost ports. The calculated flow fields and the resulting fuel distributions are shown to 
be markedly different with the different geometries. The calculated results indicate that: 
velocity profiles vary with time and are not uniform across the ports; boost port flow at high 
boost angles breaks up the toroidal vortices in the cylinder that are generated by the side ports 
and puts more fuel into the cylinder head dome and; trapping efficiency increases with increased 
boost angle. These results suggest that the computational methods developed in this work will 
be useful as a design tool for assessing the influence of engine design parameters on scavenging 
efficiencies in two-stroke engines. 

3. Gonzalez D., M.A., Borman, G.L. , and Reitz, R.D. "A Study of Diesel Cold 
Starting using both Cycle Analysis and Multidimensional Calculations," 
Accepted for publication in SAE Transactions, SAE paper 910180. 

The physical in-cylinder processes and ignition during cold starting have been studied using 
computational models, with particular attention to the influences of blowby, heat transfer during 
the compression stroke, spray development, vaporization and fuel /air mixture formation and 
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ignition. Two different modeling approaches were used. A thermodynamic zero dimensional 
cycle analysis program in which the fuel injection effects were not modeled, was used to 
determine overall and gas exchange effects. The three-dimensional KIVA-II code was used to 
determine details of the closed cycle events, with modified atomization, blowby and spray/wall 
impingement models, and a simplified model for ignition. The calculations were used to obtain an 
understanding of the cold starting process and to identify practical methods for improving cold 
starting of direct injection diesel engines. It was found that, blowby gas flow represents an 
important source of reductions for the cylinder gas temperature at lower cranking speeds, opposing 
the squish flow. Overfueling and advanced injection increase the amount of fuel evaporated. The 
spray intact core is extended at low temperatures, spray wall impingement phenomena are 
characterized by low impact velocities and the bouncing of liquid drops enhances the limited 
fuel-air mixture formation. Failure to achieve successful ignition at low initial air temperatures 
was predicted by the simple kinetics model. 

4. Reitz, R.D. and Rutland, C.J. "3-D Modeling of Diesel Engine Intake Flow 
Combustion and Emissions,” SAE Paper 911789. 

Manufacturers of heavy-duty diesel engines are facing increasingly stringent, emission standards. 
These standards have motivated new research efforts towards improving the performance of 
diesel engines. The objective of the present program is to develop a comprehensive analytical 
model of the diesel combustion process that can be used to explore the influence of design changes. 
This will enable industry to predict the effect of these changes on engine performance and 
emissions. A major benefit of the successful implementation of such models is that engine 
development time and costs would be reduced through their use. The computer model is based on 
the three-dimensional KIVA-II code, with state-of-the-art submodels for spray atomization, 
drop breakup/coalescence, multi-component fuel vaporization, spray/wall interaction, ignition 
and combustion, wall heat transfer, unbumed HC and NOx formation, and soot and radiation. 
The accuracy of the predictions is assessed by comparison with available experimental data. 
Improved combustion, wall heat transfer and spray/wall impingement submodels have been 
implemented in KIVA during the first year activity. In addition, work is in progress on a revised 
atomization model, since preliminary results show that existing atomization models are 
inaccurate under conditions of high gas temperature and pressure (e.g., turbocharged conditions). 
Finally, a methodology is being developed for modeling the intake flow process to provide more 
realistic initial conditions for engine computations. 

5. M. A. Gonzalez D. and R. D. Reitz, "Modeling Diesel Engine Spray 
Vaporization and Combustion," Proceedings of ICLASS-91, NIST, Gaithersburg, 
MD., July, 1991. 

Diesel engine in-cylinder processes have been studied using computational models with 
particular attention to spray development, vaporization, fuel/air mixture formation and 
combustion in conditions of high temperature and high pressure. A thermodynamic zero- 
dimensional cycle analysis program was used to determine initial conditions for 
multidimensional calculations. A modified version of the time-dependent, three-dimensional 
computational fluid dynamics code KIVA-II, with a detailed treatment for the spray 
calculations and a simplified model for ignition, was used to determine details of the closed 
cycle events. These calculations were used to obtain an understanding of the potential 
predictive capabilities of the models. It was found that there is a strong sensitivity of the 
spray calculations to numerical grid resolution. However, if proper grid resolution is used, the 
spray calculations were found to reproduce experimental data adequately for non-vaporizing 
sprays. However, for vaporizing sprays in high temperature engine environments the 
computations underpredicted measured gas phase (vapor) penetration results substantially. 
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This underprediction of spray penetration reduces the accuracy of combustion predictions 
greatly. The atomization drop size was found to be a key parameter influencing spray 
penetration predictions and this indicates that improved atomization models are needed for 
engine conditions. 

6. Z. W. Lian and R. D. Reitz, "Effect of Vaporization and Gas 
Compressibility on Liquid Jet Breakup," Submitted for Publication, Physics of 
Fluids, 1991. 

A linear stability analysis is presented for an evaporating jet. The development of the surface 
hydrodynamic instability is assumed to be much faster than the surface evolution due to 
evaporation. The process is then considered as quasi-steady, and the normal mode method for the 
steady basic solution is applicable as an approximation. It is found that for low speed jets 
undergoing Rayleigh breakup, jet surface evaporation is a destabilizing factor while for high 
speed atomizing jets, surface evaporation becomes stabilizing. This is due to the fact that the 
evaporation flux distributions at the troughs and crests of the waves on the surface of the liquid 
jet are different for these two cases. The effect of gas compressibility is also analyzed. For 
subsonic jets, the maximum growth rate and the corresponding wavenumber is underestimated by 
neglecting the gas compressibility since the gas pressure and gas density at the interface is higher 
than that predicted by the conventional incompressible gas theory. 


7. M. A. Gonzalez D., Z. W. Lian and R. D. Reitz, "Modeling Diesel Engine 
Spray Vaporization and Combustion," Submitted for publication for the 1992 
SAE Congress and Exposition, Detroit, MI. 

Diesel engine in-cylinder combustion processes have been studied using computational models 
with particular attention to spray development, vaporization, fuel/air mixture formation and 
combustion in conditions of high temperature and high pressure. A thermodynamic zero- 
dimensional cycle analysis program was used to determine initial conditions for the 
multidimensional calculations. A modified version of the time-dependent, three-dimensional 
computational fluid dynamics code KIVA-II, with a detailed treatment for the spray 
calculations and a simplified model for combustion, was used for the computations. These 
calculations were used to obtain an understanding of the potential predictive capabilities of the 
models. It was found that there is a strong sensitivity of the results to numerical grid resolution. 
However, with proper grid resolution, the calculations were found to reproduce experimental 
data for non-vaporizing and vaporizing sprays. However, for vaporizing sprays in high 
temperature engine environments with combustion, extremely fine grids are indicated. 
Computations made with the coarse grid sizes that are typically used underpredict measured gas 
phase (vapor) penetration results substantially. This underprediction of spray penetration 
reduces the accuracy of combustion predictions greatly. A study was made of factors that cause 
the observed sensitivity of the results to grid size in highly vaporizing sprays. The atomization 
drop size and the fuel vaporization rate were found to be key parameters that influence spray 
penetration predictions. Models for the processes that influence these parameters such as 
atomization, vapor diffusion and condensation processes are discussed. 

8. Kong, S.-C., Ayoub, N., and Reitz, R.D., "Modeling Combustion in 
Compression Ignition Homogeneous Charge Engines," Submitted for 
publication for the 1992 SAE Congress and Exposition, Detroit, MI. 


The combustion mechanism in a Compression Ignition Homogeneous Charge (CIHC) engine was 
studied. Previous experiments done on a four-stroke CIHC engine were modeled using the KIVA- 
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II code with modifications to the combustion, heat transfer, and crevice flow submodels. A 
laminar and turbulent characteristic time combustion model that has been used for spark-ignited 
engine studies was extended to allow predictions of ignition. The rate of conversion from one 
chemical species to another is modeled using a characteristic time which is the sum of a laminar 
(high temperature) chemistry time, an ignition (low temperature) chemistry time, and a 
turbulence mixing time. The ignition characteristic time was modeled using data from elementary 
initiation reactions and has the Arrhenius form. It was found to be possible to match all engine 
test cases reasonably well with one set of combustion model constants. Combustion was found to be 
controlled by chemical kinetic rates up to the time of ignition. After ignition, comparisons 
between measured cylinder pressure data and predicted pressures showed that good levels of 
agreement were not possible unless the turbulent mixing time scale was included in the combustion 
model. This important result implies that turbulent mixing, flame stretch and partial extinction 
phenomena control the rate of combustion after ignition, even in this engine that is characterized 
by homogeneous mixtures and the absence of a propagating flame. Ignition is controlled at low 
temperatures by the ignition time scale. The high temperature laminar chemistry never plays a 
role in determining the combustion rate since it is generally smaller than the other characteristic 
times. The combustion regime is classified as being between the reaction sheet and distributed 
reaction combustion. 


9. Kuo, T.-W. and Reitz, R.D., "Three-Dimensional Computations of 
Combustion in Premixed-Charge and Fuel-Injected Two-Stroke Engines, " 
Submitted for publication for the 1992 SAE Congress and Exposition, Detroit, 
MI. 

Combustion and flow were calculated in a spark-ignited two-stroke crankcase-scavenged engine 
using a laminar and turbulent characteristic time combustion submodel in the three-dimensional 
KIVA code. Both premixed-charge and fuel-injected cases were examined. A multi-cylinder 
engine simulation program was used to specify initial and boundary conditions for the 
computation of the scavenging process. A sensitivity study was conducted using the premixed- 
charge engine data. The influence of different port boundary conditions on the scavenging process 
was examined. At high delivery ratios, the results were insensitive to variations in the 
scavenging flow or residual fraction details. In this case, good agreement was obtained with the 
experimental data using an existing combustion submodel, previously validated in a four-stroke 
engine study. However, at low delivery ratios, both flow-filed and combustion-model details 
were important, and the agreement with experiment was poor using the existing combustion 
submodel, which does not account for the effect of residual gas concentration. To improve the 
agreement between modeling and experimental results, a modified combustion submodel was 
introduced that includes the effect of residual gas concentration on the laminar characteristic 
time. With the new submodel, agreement with the experiment has been improved considerably 
for all cases considered in this study. These levels of agreement between experiment and 
computations are similar to those found in previous applications of the laminar and turbulent 
characteristic- time combustion submodel to four-stroke engine combustion. Further improvement 
of the combustion submodel was made difficult by the observed coupling between the in-cylinder 
flow-field and the combustion-model details at low delivery ratios. 
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APPENDIX 3 FORTRAN SUBROUTINES 


The three FORTRAN subroutines given in this appendix have been tested for 
diesel sprays (without combustion) (subroutines atomize and walint), and for 
premixed-charge spark-ignited engine combustion (subroutine chemprn), to 
date. The subroutines need to be compiled and linked together with the 
standard KTVA subroutines. Where noted they replace existing KIVA 
subroutines. 

Subroutine atomize (replaces KIVA subroutine break) 

This subroutine is called from the main program. It computes drop breakup 
by modeling the atomization process using results from a stability analysis of 
liquid jets. The model equations and implementation details are described in 
Ref [18]. 

Subroutine walint (new subroutine to compute spray/ wall impingement) 

This subroutine computes the interaction of a drop and a wall following the 
models and equations given in Refs. [15] and [16]. It should be called in 
subroutine pmove after the call to pfind for those drops that have left the 
domain (i.e., imom=10000). 

Subroutine chemprn (replaces KIVA subroutine chem) 

This subroutine is based on the laminar and turbulent characteristic time 
’Princeton' combustion model used in Refs. [11,33] and [34]. The chemistry 
model constants are currently setup for propane fuel. References [11,33 & 34] 
should be consulted for details of model constants for other fuels. The 
subroutine also contains the Zeldo'vich NOx model as given by 
J.B. Heywood, Prog Energy Comb. Sci., Vol. 1, p. 135, 1976. 
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subroutine atomize 
include 'comd.com' 
include 'part.com' 

C +++ 

C +++ ATOMIZATION MODEL - BASED ON STABILITY OF LIQUID JETS 
C +++ SEE REITZ ATOMISATION AND SPRAY TECHNOLOGY PAPER 
C +++ 

CNST1=0.188 
CNST2=10.0 
ALPHA=0.61 
NPP = NP 
C +++ 

DO 100 N=1 ,NP 
IMOM=l4MOM(N) 

TDROP=TP(N) 

IF(IMOM.GT.99999 .OR. TDROP.GE.TCRIT) GO TO 100 
RELVEL(N)=SQRT((U(IMOM)+UTRB(N)-UP(N))“2 

1 +(V(IMOM)+VTRB(N)-VP(N))“2 

2 +(W(IMOM)+WTRB(N)-WP(N))**2) 

IF(RELVEL(N).EQ.0.0) GO TO 100 

C +++ 

SURTEN= max (STM*TDROP+STB,1.E-6) 

I4=I4P(N) 

RDROP=RADP(N) 

DROPN=PARTN(N) 

WEBER =RELVEL(N)“2*RDROP/SURTEN 
WEBERG=RO(I4)*WEBER 
WEBERL=RHOP ‘WEBER 
TB=TDROP*0.1 
ITB=INT(TB) 

FR=TB-FLOAT (ITB) 

VISCP=FR‘VISLIQ(ITB+2)+(1 ,0-FR)*VISLIQ(ITB+1 ) 

VISCP= max (VISCP.1.E-10) 

XNU =VISCP/RHOP 
REYNOL *RDROP‘RELVEL(N)/XNU 
OHN - SORT (WEBERL)/REYNOL 
TAYLOR = OHN*SQRT(WEBERG) 

FREQ *= SQRT (SURTEN/(RHOP*RDROP**3)) 

C:::::: FROM JET STABILITY DISPERSION RELATIONSHIP CURVE-FIT :::::::::: 
DENOM = (1.+0.865*WEBERG**1.67)**0.6 
WAVLNG = 9.02*(1.+0.45*OHN**0.5)*(1.+0.4*TAYLOR**0.7)/DENOM 
WAVLNG = WAVLNG'RDROP 

GROWTH = (0.34+0. 385*WEBERG**1 ,5)/((1 ,+OHN)*(1 .+1 ,4*TAYLOR**0.6)) 
GROWTH = GROWTH*FREQ 

TSHATT = 3.788*CNST2*RDROP/(GROWTH*WAVLNG) 

IF (WAVLN G.LT. RDROP/ALPHA) GO TO 40 
IF (TBREAK(N).LE.O.) GO TO 100 
TBREAK(N) = TBREAK(N) + DT 
IF (TBREAK(N).LT.TSHATT) GO TO 100 
C:::::: ENLARGE DROP ONE TIME ONLY ::::::::::::::::::::::::::::::::: 

TBREAK(N) = 0. 

RADEQ2 * (3.‘PIO2*RDROP*RDROP*RELVEL(N)/GROWTH)**0.333333 

RADEG1 = (0.75*WAVLNG‘RDROP*RDROP)**0.33333 

RADP(N) = min (RADEQ1.RADEQ2) 



PARTN(N) = DROPN‘RDROP“3/RADP(N)"3 
GOTO 100 

C:::::: BREAK-UP DROP :::::::::::::::::::::::::::::::::::::::::::::::: 

40 CONTINUE 

RADEQB - ALPHA*WAVLNG 
DTSHAT = DT/TSHATT 
C +++ 

RADP(N)=(RDROP+DTSHAT*RADEQB)/(1 ,+DTSHAT) 
PARTN(N)=DROPN*RDROP**3/RADP(N)**3 

C +++ 

IF (XDROPN(N).EQ.O.) SHEDMS(N) = 0. 

SHEDMS(N) = SHEDMS(N) - DROPN*PI403R*(RADP(N)**3-RDROP**3) 
IF (XDROPN(N).EQ.O.) XDROPN(N) = DROPN 
IF (SHEDMS(N).LT.0.03*PMINJ) GO TO 100 
C :: CREATE NEW PARCEL IF IT IS BIG ENOUGH AND IT HAS ENOUGH DROPS :: 
PARTNP = SHEDMS(N)/(PI403R‘RADEQB“3) 

IF (PARTNP.LT.XDROPN(N)) GO TO 100 
PARTN(N) = XDROPN(N) 

XDROPN(N) = 0. 

NPP = NPP + 1 
XDROPN(NPP) = 0. 

TBREAK(NPP) = 0. 

RADP(NPP) = RADEQB 
RADPP(NPP) = RADP(NPP) 

PARTN(NPP) = PARTNP 

C :::::::: CHANGE DROP VELOCITY. HERE, VI IS A VECTOR ORTHOGONAL 
C +++ TO UREL, N1 IS A UNIT VECTOR IN DIRECTION OF VI , V2 
C +++ IS A 2ND VECTOR ORTHOGONAL TO UREL AND N1, AND N2 IS 
C +++ A UNIT VECTOR IN DIRECTION OF V2 
URELX=UP(N)-(U(IMOM)+UTRB(N)) 

URELY=VP(N)-(V(IMOM)+VTRB(N)) 

URELZ=WP(N)-(W(IMOM)+WTRB(N)) 

RRELVL=URELX/(RELVEL(N)**2+1 .E-20) 

VI X=1 .0-RRELVL*URELX 
VI Y=0.0-RRELVL*URELY 
V1Z=0.0-RRELVL*URELZ 

VI Y=CVMGT(1.0,V1 Y.V1X.EQ.0.0 .AND. V1Y.EQ.0.0 .AND. V1Z.EQ.0.0) 
RV1 =1 ,0/SQRT(V1 X**2+V 1 Y**2+V 1 Z**2) 

EN1X=V1X*RV1 
EN1 Y=V1 Y*RV1 
EN1Z=V1Z‘RV1 

V2X=URELY*EN1Z-URELZ*EN1 Y 
V2Y=URELZ*EN 1 X-U RELX*EN 1 Z 
V2Z=URELX*EN1 Y-URELY*EN1X 

IF(URELX.EQ.0.0 .AND. URELY.EQ.0.0 .AND. URELZ.EQ.0.0) THEN 
V2X=0.0 
V2Y=1 .0 
V2Z=0.0 
ENDIF 

RV2=1.0/SQRT(V2X‘*2+V2Y**2+V2Z**2) 

EN2X=V2X‘RV2 
EN2Y=V2Y*RV2 
EN2Z=V2Z*RV2 
C +++ 
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C +++ TO DETERMINE VELOCITY DIRECTION, INCREMENT NORMAL TO DROP PATH. 
C +++ DELV IS MAGNITUDE OF VELOCITY CHANGE 
C +++ 

DELV-CNST1 *WAVLNG*GROWTH 
THV=PI2*FRAN(0.) 

COSTHV=COS(THV) 

SINTHV=SIN(THV) 

UP(NPP)=UP(N)+DELV*{COSTHV*EN1 X+SINTHV*EN2X) 
VP(NPP)=VP(N)+DELV*(COSTHV*EN1 Y+SINTHV*EN2Y) 
WP(NPP)=WP(N)+DELV*(COSTHV*EN1Z+SINTHV*EN2Z) 

C +++ 

TP(NPP) = TP(N) 

XP(NPP) - XP(N) 

YP(NPP) - YP(N) 

ZP(NPP) = ZP(N) 

I4P(NPP) - 14 

l4MOM(NPP) = 14 
UTRB(NPP) * UTRB(N) 

VTRB(NPP) -VTRB(N) 

WTRB(NPP) -WTRB(N) 

TURBT(NPP) - TURBT(N) 

SPDRAG(NPP) = SPDRAG(N) 

RELVEL(NPP) = RELVEL(N) 

C 

C +++ DROPLET COLLISION (O’ROURKE MODEL) APPLIED TO THE SHED DROPS 
C +++ CONTAINED IN A SWEPT VOLUME = 4 PI RDROP”2 RLVL DTSHED. 

C +++ XNCOLL IS THE PROBABLE TOTAL NUMBER OF COLLISIONS BETWEEN A 
C +++ COLLECTOR DROP WITH THE OTHER N IN THE PARCEL IN TIME DTSHED. 

C +++ XNCBAR IS THE PROBABLE NUMBER OF COLLISIONS THE COLLECTOR 
C +++ DROP EXPERIENCES DURING THE TIME INTERVAL DTSHED. 

C +++ PNOCOL IS THE PROBABILITY OF NO COLLISIONS. 

C +++ 

SUMRAD = 2.‘RADP(NPP) 

C +++ XNCOLL = 1.*PARTN(NPP)*RVOL*PI*SUMRAD"2‘RLVL*DTSHED 
C+++ XNCBAR = XNCOLL/1 . 

XNCOLL - PARTN(NPP)*(SUMRAD/RDROP)**2 
XNCBAR = XNCOLL 
PNOCOL = 0. 

C PNOCOL=EXP(-XNCBAR) 

IF (XNCBAR.LT.100.) PNOCOL=EXP(-XNCBAR) 

PNOCOL= max (PNOCOL, 1.0E-06) 

XX=FRAN(0.) 

IF(XX.LE.PNOCOL) GO TO 60 

C +++ COLLISION OCCURS, ECOAL IS THE PROBABILITY OF COALESCENCE 
FGAM=1 .3 

WECOL=RHOP*RELVEL(NPP)*RELVEL(NPP)*RADP(NPP)/SURTEN 
ECOAL= min (1 „2.4‘FGAM/(WECOL+1 .E-20)) 

YY=FRAN(0.) 

IF(YY.GT. ECOAL) GO TO 60 
C +++ CASE 1 : COALESCENCE 
ZZ=PNOCOL 
VNU=XNCBAR‘PNOCOL 
DO 20 NK=2,1000 
ZZ=ZZ+VNU 
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IF(ZZ.GE.XX) GO TO 30 
VNU=VNU*XNCBAR/FLOAT(NK) 

20 CONTINUE 
30 COLLN=FLOAT(NK-1 ) 

C +++ PROHIBIT MORE COLLISIONS THAN ARE PHYSICALLY POSSIBLE 
IF(COLLN.GE.PARTN(NPP)) COLLN=PARTN(NPP)-1 . 

RATIO = PARTN(NPP)/(PARTN(NPP)-COLLN) 

XX - RATIO*RADP(NPP)**3 
RADP(NPP)=(XX)**0.3333333333333333 
60 CONTINUE 

RADPP(NPP) * RADP(NPP) 

PARTN(NPP) = SHEDMS(N)/(PI403R*RADP(NPP)**3) 
SHEDMS(N) - 0. 

SHEDMS(NPP) =0. 


100 CONTINUE 
NP * NPP 
RETURN 
END 
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subroutine walint 


c 

c include test for the weber number of drops impinging a surface, 

c based on least square aproximation to the experimental 
c results of Wachters and Westerling. 

c October, 1989. 

c 

Include 'comd.com' 
include 'part.com' 
kpsq = kho(nho) 
if (nho.eq.O) kpsq=nzp 
KPSQ1= KPSQ - 1 
c 

if (radp(npn).le.O.) return 
WALHIT(NPN)=WALHIT(NPN)+1 .0 

::: CALCULATE THE INWARD NORMAL TO THE WALL AND DEFINE (XPOINT.YPOINT, 
::: ZPOINT AS A POINT ON THE WALL. 

I4=I4P(NPN) 

I=ITAB(I4) 

J=JTAB(I4) 

K=KTAB(I4) 

11=14+1 
12=11 +NXP 
I3=I4+NXP 
I8=I4+NXPNYP 
15=18+1 
I6=I5+NXP 
I7-I8+NXP 
WPIST = WP ISTN 
IF (K.GT.KPTOP) WPIST=0. 

::: FOR K=1 
KPED = 1 
IPED = 1 

IF (I.GE.IPED) KPED = 1 
ITOP = 0 
IBOT = 0 

IF (K. EQ. KPED. OR. K. EQ. KPTOP) IB0T=1 
IF (K.EQ.NZ.0R.K.EQ.KPSQ1) ITOP=1 
IF(K.EQ.KPED. AND.F(I1 J.NE.O.O) THEN 
IF(IBOT.EQ.1.AND.F(I1).NE.O.O) THEN 

/3\ 

3-3 — / 2 
/ 3 / 

/ / 

4 1 

XPOINT=X(l2) 

YPOINT=Y(l2) 

ZP0INT=Z(I2) 

XA=X(I2)-X(I4) 

YA=Y(I2)-Y(I4) 

ZA=Z(I2)-Z(I4) 
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c 

c 

c 

c 

c 

c 

c 

c 

c 

c 


c 

c 

c 

c 

c 

c 

c 


c 

c 

c 

c 

c 

c 

c 

c 

c 

c 


XB=X(I3)-X(I1) 

YB=Y(I3)-Y(I1) 

ZB=Z(I3)-Z(I1) 

::: FOR THE BOTTOM CORNER CELL 
ELSE IF(K.EQ.KPED) THEN 
ELSE IF(IBOT.EQ.I) THEN 

~ 3 \ /3 6 

\ / 3 
V/ 53 
3-/ 3-3 2 
// 3 / 

/ 3/ 

4 1 

XPOINT=X(l6) 

YPOINT=Y(l6) 

ZPOINT=Z(l6) 

XA=X(I6)-X(I4) 

YA=Y(I6)-Y(I4) 

ZA=Z(I6)-Z(I4) 

XB=X(I3)-X(I5) 

YB=Y(I3)-Y(I5) 

ZB=Z(I3)-Z(I5) 

::: FOR TOP WALL K=NZ OR K=KPSQ1 
ELSE IF(ITOP.EQ.1.AND.F(I1).NE.O.O) THEN 

7 / 6 

/ 3 / 

/ / 

8 3 5 

\3/ 

XPOINT=X(l6) 

YP0INT=Y(I6) 

ZP0INT=Z(I6) 

XA=X(I6)-X(I8) 

YA=Y(I6)-Y(I8) 

ZA=Z(I6)-Z(I8) 

XB=X(I5)-X(I7) 

YB=Y(I5)-Y(I7) 

ZB=Z(I5)-Z(I7) 

::: FOR THE TOP CORNER CELL 
ELSE IF(ITOP.EQ.I) THEN 

7 

/ /3 6 

/ / 3 

8 / / 53 

\ 3—3 2 

/\ 3 / 

/ 3/ 

3/ 3/1 


XPOINT=X(l2) 

YPOINT=Y(l2) 

ZPOINT-Z(12) 
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XA=X(I2)-X(I8) 

YA-Y(I2)-Y(I8) 

ZA=Z(I2)-Z(I8) 

XB=X(I1)-X(I7) 

YB*Y(I1)-Y(I7) 

ZB=Z(I1)-Z(I7) 

C ::: FOR CONSTANT I WALL 
ELSE 
C 

C /3 6 

C 5/3 

C < 3 - 3 

C 3 — 3-3 2 

C / 3/ 

C / 3/ 

C 4 1 

XPOINT=X(l5) 

YPOINT=Y(l5) 

ZP0INT=Z(I5) 

XA=X(I5)-X(I2) 

YA=Y(I5)-Y(I2) 

ZA=Z(I5)-Z(I2) 

XB=X(I6)-X(I1) 

YB=Y(I6)-Y(I1) 

ZB=Z(I6)-Z(I1) 

END IF 

XNORM=YA*ZB-ZA‘YB 
YNORM=ZA*XB-XA‘ZB 
ZNORM=XA*YB-YA*XB 

SNORM=SQRT(XNORM**2+YNORM**2+ZNORM**2) 

XNORM=XNORM/SNORM 
YNORM=YNORM/SNORM 
ZNORM=ZNORM/SNORM 

UPREL=UP(NPN) 

VPREL=VP(NPN) 

WPREL=WP(NPN)-WPIST 

QDOTN=UPREL*XNORM+VPREL*YNORM+WPREL*ZNORM 

::: GAMMAA IS THE ANGLE BETWEEN THE NORMAL AND RELATIVE VELOCITY 

QPARTs=SQRT(UPREL**2+VPREL**2+WPREL**2) 
GAMMAA=acos(-(UPREL*XNORM+VPREL*YNORM+WPREL*ZNORM)/QPART) 

C ::: DETERMINE THE POINT (XWALHT.YWALHT.ZWALHT) THE DROP HIT THE WALL. 
C ::: DPOINT IS THE NORMAL DISTANCE FROM THE PRESENT DROP TO THE WALL. 

C ::: 

DPOINT=(XP(NPN)-XPOINT)*XNORM+(YP(NPN)-YPOINT)*YNORM+ 

1 (ZP(NPN)-ZPO!NT)*ZNORM 
DPOQCG=DPOINT/(QPART*cos(GAMMAA)) 

XWALHT =XP(NPN)+DPOQCG*UPREL 
YWALHT *YP(NPN)+DPOQCG*VPREL 
ZWALHT =ZP(NPN)+DPOQCG*WPREL 
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surten=max(stm‘tp(npn)+stb,1 ,0e-06) 
wedrop=rhop*(qdotn**2)*2.0*radp{npn)/surten 
c 

c test to determine the nature of the impact 
c 

if(wedrop.le.80.0) then 
c weber number after impact 

wedropo=0.67852227*wedrop*(exp((-4.41 5135e-02)*wedrop)) 
c normal velocity after impact 

qdotn2=qdotn*(1 .0+sqrt(wedropo/wedrop)) 
c velocity components 

up(npn)=uprel-qdotn2*xnorm 

vp(npn)=vprel-qdotn2*ynorm 

wp(npn)=wprel-qdotn2*znorm+wpist 

c 

c new location 

xp(npn)=xwa!ht+abs(dpoqcg)*up(npn) 
yp(npn)=ywalht+abs(dpoqcg)‘vp(npn) 
zp(npn)=zwalht+abs(dpoqcg)*wp(npn) 

distwp=sqrt((xp(npn)-xwalht)*‘2 + (yp(npn)-ywalht)**2+ 

* (zp(npn)-zwalht)**2) 

else 

::: THE TANGENTIAL VECTOR ON THE WALL SURFACE IN THE DIRECTION OF 
::: THE RELATIVE VELOCITY VECTOR IS GIVE BY (XTAN,YTAN,ZTAN) 

XT AN=UPREL-QDOTN‘XNORM 
YT AN=VPREL-QDOTN* YNORM 
ZTAN=WPREL-QDOTN*ZNORM 

STAN=SQRT(XTAN**2+YTAN**2+ZTAN"2) 

XTAN=XT AN/STAN 
YT AN = YT AN/ST AN 
ZTAN=ZTAN/STAN 

::: THE BINORMAL VECTOR B=(T X N) IS IN THE PLANE OF THE SURFACE 
::: B=(XBIN,YBIN,ZBIN) 

XBIN=YTAN*ZNORM-ZTAN*YNORM 
YBIN=ZTAN*XNORM-XTAN*ZNORM 
ZBIN=XTAN*YNORM-YTAN*XNORM 

C ::: 

C ::: DETERMINE PSI THE ANGLE RELATIVE TO THE TANGENTIAL VECTOR IN 
C ::: THE PLANE OF THE WALL TO MOVE THE DROP. BETA IS THE PARAMETER 
C ::: DETERMINED FROM THE ANGLE GAMMAA. 

C ::: 

BETA=PI‘SQRT(sin(GAMMAA)/(1.0-sin(GAMMAA))) 

XXX=FRAN(0.0) 

YYY=FRAN(0.0) 

PSI=-(PI/BETA)* log(1 .0-XXX*(1 .0-EXP(-BETA))) 

IF(YYY.GT.0.5) PSI=-PSI 
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::: FIND THE DIRECTION OF THE NEW RELATIVE VELOCITY 

COSPSI=cos(PSI) 

SINPSI=sin(PSI) 

XVNEW=XTAN*COSPSI+XBIN*SINPSI 
YVNEW=YTAN*COSPSI+YBIN*SINPSI 
ZVNEW=ZTAN*COSPSI+ZBIN*SINPSI 

::: SET THE NEW RELATIVE VELOCITY AS SOME FRACTION (FRACT) OF THE OLD 
::: RELATIVE VELOCITY AND UP DATE THE NEW ABSOLUTE VELOCITIES. 

FRACT=1 .0 

VELNEW=FRACT*QPART 
UP(NPN)=VELNEW*XVNEW 
VP(NPN)=VELNEW*YVNEW 
WP(NPN)=VELNEW*ZVNEW + WPIST 
c 

c new location 

DMOVE=2.0*RADP(NPN) 

XP(NPN)=XWALHT+DMOVE*XNORM 

YP(NPN)=YWALHT+DMOVE*YNORM 

ZP(NPN)=ZWALHT+DMOVE*ZNORM 

c 

END IF 
C ::: 

RETURN 

END 
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subroutine chemprn 
include 'comd.com' 

C +++ 

CROLF.1/88... PRINCETON COMBUSTION MODEL 

C ... FROM ABRAHAM ET AL. COMBUSTION & FLAME, VOL 60, P. 309, 1985. 

C +++ 

c DIMENSION YY(9),CK(9),A(5),B(5),D(5) 

DIMENSION YY(12),CK(12),A(5),B(5),D(5) 

REAL KC02.KH20 

F0(Q)=1 6.31 8*0-18.972 

FI (Q)=1 1 .221 3*Q*‘3-21 .8752*0**2+1 4.9469*0-2.654 
F2(Q)=0.9952*Q-0.351 4 
DATA EPSIT.ITMAX/1 .E-03,90/ 

DATA BCO2,ECO2,SCO2/5.86E-06, 68353.9, 0.9908/ 

DATA BH20,EH20,SH20/8.60534E-06, 1041 5.9,0.9894/ 

DATA IFUEL,I02,IN2,IC02,IH20,IH2,IC0,IN0/ 

41.2,3, 4, 5, 7, 11, 12/ 

DATA CM1, CM2 , CM3 / 

& 8.25,0.142,0.235/ 

C NOTE: CM1 INCREASED TO FOR PR=0.66. USE 1 .9 FOR SHORT SPARK ELSE 1 .8 
data taud,hspark/0.,0.3/ 

DATA ANUMC,ANUMH,ANUMO/3.,8.,0./ 

data BF,EA,SSL1 ,SSL2,SSL3,ALPHA1 ,ALPHA2,BETA1 .BETA2.ACTEX1 , ACTEX2/ 
& 1 ,54e-12,3.0e04, 34.22,1 38.65,1 .08,2.1 8, 0.8, 0.1 6,0.22,0.08,1 . 1 5/ 

data temspk,pspark,gamspk,resid,phispk/0.,0.,0.,0.,0./ 

TCHEM = 1.0E-10 
nspark=2 

if (iignl(2).eq.O) nspark=1 
C %%% LOCATE THE IGNITION CELL 
IF (NHO.GT.O) GO TO 191 
DO 190 N=1,NSPARK 
DO 180 K=1,NZ 
I4K=(K-1 )*NXPNYP 
WJK=I4K+(JIGNF(N)-1 )*NXP 
I4=I4JK+IIGNL(N) 

IF(ZHEAD-Z(I4).LE.HSPARK) GO TO 185 
180 CONTINUE 
185 KSPARK = K 

IF (KSPARK.GE.NZ) KSPARK = NZ-1 

KIGNB(N)=KSPARK 

KIGNT (N)=KSPARK 

190 CONTINUE 

191 CONTINUE 

I4R = (KIGNB(1 )-1 )*NXPNYP + (JIGNF(1)-1)*NXP + IIGNL(I) 

COMB STORE UNBURNED GAS PARAMETERS FOR IGNITION MODEL 

IF (TAUD.NE.0.) GO TO 1000 
PHISPK=0. 

DO 900 K=1,NZ 
I4B=(K-1 )‘NXPNYP 
DO 800 J=1 ,NY 
I4=I4B+(J-1)*NXP+1 
DO 700 1=1, NX 
IF(F(I4).EQ.0.) GO TO 700 
DO 200 ISP=1 ,nsp 
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YY(1SP)=SPD(I4, ISP)/RO(l4) 

CK(ISP)=SPD(I4,ISP)*RMW(ISP)/R0(I4) 

200 CONTINUE 

0XYGN=2.*(CK(I02)+CK(IC02))+CK(IH20)+CK(IC0)+ANUM0*CK(IFUEL) 
HYDRN=2.*(CK(IH20)+CK(IH2)) +ANUMH*CK(IFUEL) 

CARBN=CK(IC02)+CK(ICO) +ANUMC*CK(IFUEL) 

IF (I4.EQ.I4R) PHISPK=CK(IFUEL)*(2.*ANUMC+0.5*ANUMH)/ 

& (ANUM0*CK(IFUEL)+2.*CK(I02)) 

SUMY=0. 

DO 210 ISP = 4,11 
210 SUMY - SUMY + YY(ISP) 
if (i4.eq.l4r) resid=sumy 

C +++ SAVE RESIDUAL CONCENTRATIONS AT THE TIME OF SPARK IN SPD(NSP) 
c SPD(I4,NSP)= SUMY*R0(I4) 

700 14=14+1 
800 CONTINUE 
900 CONTINUE 
C +++ 

TEMSPK = TEMP(I4R) 

PSP ARK = P(I4R) 

GAMSPK = GAMMA(I4R) 

C 

1000 CONTINUE 

ESTIMATE THE UNBURNED GAS TEMPERATURE 

AND LAMINAR FLAME SPEED NEAR THE SPARK CELL 

TEMPI4 = TEMSPK*(P(I4R)/PSPARK)**(1.-1. /GAMSPK) 

SSUBL0 = SSL1 - SSL2*(PHISPK-SSL3)**2 
IF (SSUBL0.LT.1.0) SSUBL0 = 1.0 
ALPHA = ALPHA1 - ALPHA2*(PHISPK - 1.) 

BETA =-BETA1 + BETA2*(PHISPK - 1.) 

SSUBL - SSUBL0*((TEMPI4/298.)** ALPHA) 
&*((P(I4R)/1.013E06)**BETA)*(1.-2.1*RESID) 

IF (SSUBL.LT.1.0) SSUBL = 1.0 
ELT = 0.1 6431 68‘(TKE(I4R)**1.5)/EPS(I4R) 

TAUD= CM1*ELT/SSUBL 

END SETUP FOR FLAME KERNEL GROWTH MONITORING 

DO 90 K=1,NZ 
I4B=(K-1 )*NXPNYP 
DO 80 J=1,NY 
I4=I4B+(J-1 )*NXP+1 
DO 70 1=1, NX 
IF(F(I4).EQ.0.) GO TO 70 
C +++ 

C +++ CONCENTRATIONS CK, MASS FRACTIONS YY 
C +++ 

DO 20 1SP=1 ,1 1 
YY(ISP)=SPD(l4,ISP)/RO(l4) 

CK(ISP)=SPD(l4,ISP)*RMW(ISP)/RO(l4) 

20 CONTINUE 
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0XYGN=2.*(CK(I02)+CK(IC02))+CK(IH20)+CK(IC0)+ANUM0‘CK(IFUEL) 
HYDRN=2.*(CK(IH20)+CK(IH2)) +ANUMH*CK(IFUEL) 

CARBN=CK(IC02)+CK(IC0) +ANUMC‘CK(IFUEL) 

EQPHI-(2.*CARBN+0.5‘HYDRN)/OXYGN 

C APPROX. TREATMENT FOR PHI>3.0 AND PHI<40 

BOT=2.*ANUMC+0.5*ANUMH-3.*ANUMO 
FUEL* 0. 

IF (EQPHI.GE.3.0) FUEL=(2.*CARBN+0.5‘HYDRN-3.*OXYGN)/BOT 
CARBN=CARBN-ANUMC*FUEL 
HYDRN=HYDRN-ANUMH*FUEL 
OXYGN=OXYGN-ANUMO*FUEL 

.... ITERATE TO FIND EQUILIBRIUM SPECIES DENSITIES 

.... USE C0C02 RATIO FROM NASA CODE FIT AS INITIAL GUESS ... 

IPHI=EQPHI 

IF (IPHI.EQ.O) COCO2=F0(EQPHI) 

IF (IPHI.EQ.1) COC02=F1 (EQPHI-1 .) 

IF (IPHI.GT.1) C0C02=F2(EQPHI) 

C0C02=2.‘EXP(C0C02) 

COCO0=COCO2 

5 EAA=1 .987*TEMP(I4) 

C20=0. 

ECEA=EC02/EAA 
IF (ECEA.GT.50.) GO TO 6 

C ... KC02 = EQUILIB. CONSTANT FOR CO + .502 = C02 

KC02=BC02*(TEMP(I4)**SC02)*EXP(ECEA) 

C20=2./(KCO2*KCO2*RO(l4)) 

C ... KH20 = KC02 / KH20' (KH20* EQ. CONST. FOR H2 + 0.5O2 = H20 ) .. 

6 KH20=BH20*(TEMP(I4)**SH20)*EXP(EH20/EAA) 
A(5)=KH20*(CARBN-0XYGN) 

A(4)=CARBN‘(2.‘KH20+1 .)+0.5*HYDRN-OXYGN*(KH2O+1 .) 
A(3)=2.*CARBN+0.5*HYDRN-OXYGN+KH2O*C20 
A(2)=(KH2O+1.)*C20 
A(1 )=C20 
B(5)=A(5) 

D(5)=B(5) 

C .. SOLVE QUARTIC EQUATION FOR CO-C02 RATIO WITH NEWTON POLY. SOLVER 
DO 8 ICONC=1,ITMAX 
DO 7 LL=1 ,4 
L=5-LL 

B(L)=A(L)+COC02*B(L+1 ) 

7 D(L)=B(L)+C0C02*D(L+1 ) 

COC21 =COC02-B(1 )/D(2) 

IF (ABS(COC21-COCO2).LT.EPSIT*COCO0) GO TO 9 
COC02=COC21 

8 IF (COCO2.LT.0.) C0C02=-C0C02 
ERR0=0. 

WRITE(6,14) ITMAX,TEMP(I4) 

14 FORMAT(1H ,2X,'C0-C02 ITERATION FAILED AFTER’,15,' ITERATIONS’, 

$’ TEMPERATURE', G1 0.4) 

COCO2=COCO0 

C ... COMPUTE LOCAL EQUILIBRIUM COMPOSITION, YY’S 

9 H2O=0.5*HYDRN/(1 .+KH20*C0C02) 
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C02=CARBN/(1 .+C0C02) 

YY(IH20)=H20*MW(IH20) 

YY(IC02)=C02*MW(IC02) 

YY(IC0)=C0C02*C02*MW(IC0) 

YY(IH2)=H20‘KH20*C0C02*MW(IH2) 

YY(I02)=1.-(YY(IH20)+YY(IC02)+YY(IC0)+YY(IH2)+YY(IN2)) 

IF (EQPHI.GE.2.) YY(IO2)=0.5*C20‘MW(IO2)/(COCO2*COCO2) 
YY(IFUEL)=1.-(YY(I02)+YY(IH20)+YY(IC02)+YY(IC0)+YY(IH2)+YY(IN2)) 
c 

DELTS = (T-T1IGN)/TAUD 
FACSPK * 1 .-EXP(-DELTS) 

C ... FIND CHARACTERISTIC REACTION TIME, TAUCHM 

CM23 - CM2/CM3 
DELYPS = 0. 

DO 55 ISP = 4,11 

55 DELYPS - DELYPS + SPD(I4,ISP)/R0(I4) 
cresIdDELYPS - DELYPS - SPD(l4,NSP)/RO(l4) 
c dypsav=delyps 

DELYPS = DELYPS - resld 

DELYFO = (SPD(I4,IFUEL)+SPD(I4,I02))/R0(I4)-(YY(IFUEL)+YY(I02)) 
DELYPS = ABS(DELYPS) 

DELYFO * ABS(DELYFO) 

FACPHI - 1.0E-20 

IF (DELYFO.NE.O..AND.DELYPS.NE.O.) FACPHI = CM23‘DELYPS/DELYFO 
TAUTRB = CM2*TKE(I4)/EPS(I4) 

CROLF IF (FACPHI.LT.1..AND.DELYFO.NE.O.) TAUTRB = TAUTRB/FACPH I 
IF (FACPHI.LT.1 ..AND.DELYFO.NE.O..AND.DELYPS.NE.O.) 

1 TAUTRB = TAUTRB/FACPHI 

C ... LAMINAR FLAME KINETIC TIME . 

EQPHI=CK(IFUEL)*(2.*ANUMC+0.5*ANUMH)/ 

& (ANUM0*CK(IFUEL)+2.*CK(I02)) 

BETA --BETA1 + BETA2*(EQPHI - 1 .) 

PREEXP - BF*TEMP(I4)*((1 .01 3E06/P(I4))“(1 .+2/BETA)) 

CKUO PREEXP = BF*TEMP(I4)*((1 .01 3E06/P(I4))‘*0.75) 

ACTEXP - (1 .+ACTEX1 *ABS(EQPHI-ACTEX2))*EA/(1 .987*TEMP(I4)) 
taulam * 0. 

if (abs(actexp).le.200.) 

&TAULAM = PREEXP*EXP(ACTEXP) 

TAULAM « 0.62‘TAULAM/(1 .-2.1 *RESID)**2 
c 

TAUCHM = TAULAM + TAUTRB'FACSPK 

C FIRST ORDER ACCURATE STIFF SCHEME 

FACEXP . 1 .-EXP(-DT/TAUCHM) 

DO 60 ISP=1,11 

RRATE - -(SPD(I4,ISP) - RO(l4)*YY(ISP))*FACEXP 

SPD(I4,ISP) = SPD(I4,ISP) + RRATE 

DECHEM = - RRATE*HTFORM(ISP)*RMW(ISP)/RO(l4) 

DECHK - ABS(DECHEM/SIE(I4)) 

SIE(I4) - SIE(I4) + DECHEM 
TCHEM = max (TCHEM, DECHK) 

C END COMBUSTION MODEL 

60 CONTINUE 
C 
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EXTENDED ZELDOVICH NOX MODEL 

... FROM J.B. HEYWOOD, PROG. ENERGY COMB. SCI., VOL. 1, P. 135, 1976 

NO = 1 

IF (NO.EQ.O) GO TO 70 
TA=TEMP(I4)*0.001 
TALOG= log(TA) 

C EQUILIBRIUM CONSTANTS FROM M&M GMR-4361 : 

OKI = 0.794709* log(TA)-1 13.208/TA+3.1 6837-0.4438 14*TA 
1 +2.69699E-02*TA**2 

QK2 - 0.431310* log(TA)- 59.655/TA+3.50335-0.340016*TA 
1 +1 .5871 5E-02*TA**2 

QK4 * 0.990207* log(TA)- 51 ,792/TA+0.99307-0.343428*TA 
1 +1 .1 1668E-02*TA**2 

QK5 - 0.652939* log(TA)+ 9.823/TA-3.93033-0.163490*TA 
1 +1.42865E-02*TA“2 

OKI * EXP(QKI) 

QK2 = EXP(QK2) 

QK4 = EXP(QK4) 

QK5 = EXP(QK5) 

C NOW DO NO KINETICS 

OXYGEN = ABS(SPD(I4,I02))*RMW(I02) 

XNITRO = ABS(SPD(I4,IN2))*RMW(IN2) 

HYDROG = ABS(SPD(I4,IH2))*RMW(IH2) 

QKNO = 4.4742*EXP(-75.677/TA) 

QK02N2 = SORT (OKI *QK2*OXYGEN*XN ITRO) 

ALFAP = QKN0*RMW(IN0)/QK02N2 
ALFA = SPD(l4,INO)‘ALFAP 
SQK1N2 = SQRT(QK1 ‘XNITRO) 

RR1 - 1.626E13*QK02N2*SQK1N2/QKN0 

RR2 = 6.625E09*TEMP(I4)*EXP(-3150./TEMP(I4))*SQK1N2*OXYGEN 
RR3 = 4.21 6E1 3*SQRT(QK1 *XNITRO*OXYGEN*HYDROG/QK5) 

FORWD - DT*2.*MW(INO)*RR1/(1.+ALFA*RR1/(RR2+RR3)) 

TOP = SPD(l4,INO)+FORWD 
BOT = O.5*(1.+SQRT(1.+4.‘TOP‘FORWD*ALFAP**2)) 

SPD(l4,INO) = TOP/BOT 

CROLF END NOX MODEL 

70 14=14+1 
80 CONTINUE 
90 CONTINUE 
RETURN 
END 
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